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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
|
/**
\file error_estimator.hpp
\brief define the function used to estimate the error boundaries of a fit
\author Junhua Gu
*/
#ifndef ERROR_EST
#define ERROR_EST
#include <core/fitter.hpp>
#include <core/freeze_param.hpp>
#include <interface/type_depository.hpp>
#include <math/num_diff.hpp>
#include <iostream>
#include <vector>
#include <string>
namespace opt_utilities
{
template <typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
typename element_type_trait<Tp>::element_type estimate_error_hessian(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,const Ts& dchi)
{
size_t porder=fit.get_param_order(pname);
holder<param_modifier<Ty,Tx,Tp,Tstr> > h;
try
{
h.reset(fit.get_param_modifier().clone());
}
catch(const param_modifier_not_defined&)
{
h.reset(0);
}
fit.clear_param_modifier();
Tp p=fit.get_all_params();
Ts a=hessian(fit.get_statistic(),p,porder,porder)/2;
if(h.get())
{
fit.set_param_modifier(*(h.get()));
}
fit.fit();
return std::sqrt(dchi/a);
}
/**
\brief calculate the error boundary of a fit, according to the given delta statistic.
\param fit the fitter that has a sucessful fit result
\param pname the name of the parameter, the error of which will be estimated
\param lower input as the initial value of the lower boundary, and output as the final result of the lower boundary
\param upper input as the initial value of the upper boundary, and output as the final result of the upper boundary
\param dchi the delta statistic corresponding to a certain confidence level
\param precision determine how precise the error bounds should be determined
*/
template <typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
void estimate_error_directly(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,typename element_type_trait<Tp>::element_type& lower,typename element_type_trait<Tp>::element_type& upper,const Ts& dchi,const Ts& precision)
{
typedef typename element_type_trait<Tp>::element_type Tpe;
std::vector<Tstr> pnames;
std::vector<Tpe> pvalues;
//Make sure we start from an optimal parameter set
fit.fit();
//stores origin parameter values
for(size_t i=0;(size_t)i<fit.get_num_params();++i)
{
pnames.push_back(fit.get_param_info(i).get_name());
pvalues.push_back(fit.get_param_info(i).get_value());
}
//ensure that the interesting parameter is free
if(fit.report_param_status(pname)=="frozen")
{
return;
}
try
{
freeze_param<Ty,Tx,Tp,Tstr>* pfp=dynamic_cast<freeze_param<Ty,Tx,Tp,Tstr>*>(&fit.get_param_modifier());
if(pfp==0)
{
assert(0);
}
*pfp+=freeze_param<Ty,Tx,Tp,Tstr>(pname);
}
catch(const param_modifier_not_defined&)
{
fit.set_param_modifier(freeze_param<Ty,Tx,Tp,Tstr>(pname));
}
//get current statistic value
Tpe current_value=
fit.get_param_value(pname);
//initial lower boundary should be a worse parameter,
//so do the upper boundary
if(lower>=current_value||upper<=current_value)
{
std::cerr<<"Error, initial lower and upper limits should be smaller and larger than the best fit value, respectively"<<std::endl;
return;
}
Ts current_chi=fit.get_statistic_value();
Tpe upper1=current_value;
Tpe step=upper-current_value;
//find a new upper bound, whose statistic is worse than the target statistic
while(1)
{
fit.set_param_value(pname,upper);
fit.fit();
Ts chi=fit.get_statistic_value();
std::cerr<<upper<<"\t"<<chi<<std::endl;
if(chi>current_chi+dchi)
{
break;
}
upper+=step;
}
//perform the bi-search
Ts target_chi=current_chi+dchi;
while(upper-upper1>std::abs(precision))
{
fit.set_param_value(pname,upper1);
fit.fit();
Ts s1=fit.get_statistic_value();
fit.set_param_value(pname,upper);
fit.fit();
Ts s=fit.get_statistic_value();
fit.set_param_value(pname,(upper+upper1)/2);
fit.fit();
Ts s01=fit.get_statistic_value();
std::cerr<<s1-target_chi<<"\t"<<s01-target_chi<<"\t"<<s-target_chi<<std::endl;
if((s1-target_chi)*(s01-target_chi)<0&&(s-target_chi)*(s01-target_chi)>=0)
{
upper=(upper+upper1)/2;
}
else if((s1-target_chi)*(s01-target_chi)>0&&(s-target_chi)*(s01-target_chi)<=0)
{
upper1=(upper+upper1)/2;
}
else
{
std::cerr<<"strange statistic structure"<<std::endl;
break;
}
}
for(size_t i=0;i<fit.get_num_params();++i)
{
fit.set_param_value(pnames[i],pvalues[i]);
}
fit.fit();
////
Tpe lower1=current_value;
step=current_value-lower;
while(1)
{
fit.set_param_value(pname,lower);
fit.fit();
Ts chi=fit.get_statistic_value();
std::cerr<<lower<<"\t"<<chi<<std::endl;
if(chi>current_chi+dchi)
{
break;
}
lower-=step;
}
fit.set_param_value(pname,current_value);
fit.fit();
while(lower1-lower>std::abs(precision))
{
fit.set_param_value(pname,lower1);
fit.fit();
Ts s1=fit.get_statistic_value();
fit.set_param_value(pname,lower);
fit.fit();
Ts s=fit.get_statistic_value();
fit.set_param_value(pname,(lower+lower1)/2);
fit.fit();
Ts s01=fit.get_statistic_value();
std::cerr<<s1-target_chi<<"\t"<<s01-target_chi<<"\t"<<s-target_chi<<std::endl;
if((s1-target_chi)*(s01-target_chi)<0&&(s-target_chi)*(s01-target_chi)>=0)
{
lower=(lower+lower1)/2;
}
else if((s1-target_chi)*(s01-target_chi)>0&&(s-target_chi)*(s01-target_chi)<=0)
{
lower1=(lower+lower1)/2;
}
else
{
std::cerr<<"strange statistic structure"<<std::endl;
break;
}
}
//restore the param_modifier
dynamic_cast<freeze_param<Ty,Tx,Tp,Tstr>& >(fit.get_param_modifier())-=freeze_param<Ty,Tx,Tp,Tstr>(pname);
//restore the origin param values
for(size_t i=0;i<fit.get_num_params();++i)
{
fit.set_param_value(pnames[i],pvalues[i]);
}
fit.fit();
}
template <typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
void estimate_error(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,typename element_type_trait<Tp>::element_type& lower,typename element_type_trait<Tp>::element_type& upper,const Ts& dchi,const Ts& precision)
{
typename element_type_trait<Tp>::element_type e=estimate_error_hessian(fit,pname,dchi);
lower=fit.get_param_value(pname)-e*2;
upper=fit.get_param_value(pname)+e*2;
estimate_error_directly(fit,pname,lower,upper,dchi,precision);
}
}
#endif
//EOF
|
