aboutsummaryrefslogtreecommitdiffstats
path: root/statistics/cstat.hpp
blob: f1c1d8f9bdb888c6b8af101aa09878df132c80c6 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
/**
   \file cstat.hpp
   \brief maximum-liklihood statistic
   \author Junhua Gu
 */

#ifndef CSTAT_HPP
#define CSTAT_HPP
#define OPT_HEADER
#include <core/fitter.hpp>
#include <math/vector_operation.hpp>
#include <iostream>
#include <cassert>


using std::cout;using std::endl;
namespace opt_utilities
{

  /**
     \brief c-statistic, max-likelihood method
     \tparam Ty the return type of model
     \tparam Tx the type of the self-var
     \tparam Tp the type of model parameter
     \tparam Ts the type of the statistic
     \tparam Tstr the type of the string used
   */
  template<typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
  class cstat
    :public statistic<Ty,Tx,Tp,Ts,Tstr>
  {
  private:
    bool verb;
    int n;
  public:
    cstat()
      :verb(true)
    {}

    void verbose(bool v)
    {
      verb=v;
    }

    const char* do_get_type_name()const
    {
      return "maximum likelihood";
    }

  public:

    statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const
    {
      // return const_cast<statistic<Ty,Tx,Tp>*>(this);
      return new cstat<Ty,Tx,Tp,Ts,Tstr>(*this);
    }

    Ts do_eval(const Tp& p)
    {
      Ts result(0);
      for(int i=(this->get_data_set()).size()-1;i>=0;--i)
	{
	  Ty model_y=this->eval_model(this->get_data_set().get_data(i).get_x(),p);
	  result-=contract(this->get_data_set().get_data(i).get_y(),std::log(model_y),result);
	}

      if(verb)
	{
	  n++;
	  if(n%10==0)
	    {
	      cout<<"a:"<<result<<"\t";
	      for(size_t i=0;i<get_size(p);++i)
		{
		  cout<<get_element(p,i)<<",";
		}
	      cout<<endl;
	    }
	  
	}
      return result;
    }
  };

    /**
     \brief c-statistic, max-likelihood method
     \tparam Ty the return type of model
     \tparam Tx the type of the self-var
     \tparam Tp the type of model parameter
     \tparam Ts the type of the statistic
     \tparam Tstr the type of the string used
   */
  template<typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
  class cstat1
    :public statistic<Ty,Tx,Tp,Ts,Tstr>
  {
  private:
    bool verb;
    int n;
  public:
    cstat1()
      :verb(true)
    {}

    void verbose(bool v)
    {
      verb=v;
    }

    const char* do_get_type_name()const
    {
      return "maximum likelihood";
    }

  public:

    statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const
    {
      // return const_cast<statistic<Ty,Tx,Tp>*>(this);
      return new cstat1<Ty,Tx,Tp,Ts,Tstr>(*this);
    }

    Ts do_eval(const Tp& p)
    {
      Ts result(0);
      for(int i=(this->get_data_set()).size()-1;i>=0;--i)
	{
	  Ty model_y=this->eval_model(this->get_data_set().get_data(i).get_x(),p);
	  result-=contract1(this->get_data_set().get_data(i).get_y(),std::log(model_y),result);
	}

      if(verb)
	{
	  n++;
	  if(n%10==0)
	    {
	      cout<<"a:"<<result<<"\t";
	      for(size_t i=0;i<get_size(p);++i)
		{
		  cout<<get_element(p,i)<<",";
		}
	      cout<<endl;
	    }
	  
	}
      return result;
    }
  };
}

#endif
//EOF