ROOT   Reference Guide
FumiliStandardMaximumLikelihoodFCN.cxx
Go to the documentation of this file.
1// @(#)root/minuit2:$Id$
2// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
10
12
13#include <vector>
14#include <cmath>
15#include <float.h>
16
17namespace ROOT {
18
19 namespace Minuit2 {
20
21
22//#include <iostream>
23
24std::vector<double> FumiliStandardMaximumLikelihoodFCN::Elements(const std::vector<double>& par) const {
25
26 // calculate likelihood element f(i) = pdf(x(i))
27 std::vector<double> result;
28 double tmp1 = 0.0;
29 unsigned int fPositionsSize = fPositions.size();
30
31
32 for(unsigned int i=0; i < fPositionsSize; i++) {
33
34 const std::vector<double> & currentPosition = fPositions[i];
35
36 // The commented line is the object-oriented way to do it
37 // but it is faster to do a single function call...
38 //(*(this->getModelFunction())).SetParameters(par);
39 tmp1 = (*(this->ModelFunction()))(par, currentPosition);
40
41 // std::cout << " i = " << i << " " << currentPosition[0] << " " << tmp1 << std::endl;
42
43 result.push_back(tmp1);
44
45 }
46
47
48 return result;
49
50}
51
52
53
54const std::vector<double> & FumiliStandardMaximumLikelihoodFCN::GetMeasurement(int index) const {
55 // Return x(i).
56 return fPositions[index];
57
58}
59
60
62 // return size of positions (coordinates).
63 return fPositions.size();
64
65}
66
67
68void FumiliStandardMaximumLikelihoodFCN::EvaluateAll( const std::vector<double> & par) {
69 // Evaluate in one loop likelihood value, gradient and hessian
70
71 static double minDouble = 8.0*DBL_MIN;
72 static double minDouble2 = std::sqrt(8.0*DBL_MIN);
73 static double maxDouble2 = 1.0/minDouble2;
74 // loop on the measurements
75
76 int nmeas = GetNumberOfMeasurements();
78 std::vector<double> & h = Hessian();
79 int npar = par.size();
80 double logLikelihood = 0;
82 h.resize( static_cast<unsigned int>(0.5 * npar* (npar + 1) ) );
84 h.assign(static_cast<unsigned int>(0.5 * npar* (npar + 1) ) , 0.0);
85
86 const ParametricFunction & modelFunc = *ModelFunction();
87
88 for (int i = 0; i < nmeas; ++i) {
89
90 // work for one-dimensional points
91 const std::vector<double> & currentPosition = fPositions[i];
92 modelFunc.SetParameters( currentPosition );
93 double fval = modelFunc(par);
94 if (fval < minDouble ) fval = minDouble; // to avoid getting infinity and nan's
95 logLikelihood -= std::log( fval);
96 double invFval = 1.0/fval;
97 // this method should return a reference
99
100 // calc derivatives
101
102 for (int j = 0; j < npar; ++j) {
103 if ( std::fabs(mfg[j]) < minDouble ) {
104 // std::cout << "SMALL values: grad = " << mfg[j] << " " << minDouble << " f(x) = " << fval
105 // << " params " << j << " p0 = " << par[0] << " p1 = " << par[1] << std::endl;
106 if (mfg[j] < 0)
107 mfg[j] = -minDouble;
108 else
109 mfg[j] = minDouble;
110 }
111
112 double dfj = invFval * mfg[j];
113 // to avoid summing infinite and nan later when calculating the Hessian
114 if ( std::fabs(dfj) > maxDouble2 ) {
115 if (dfj > 0)
116 dfj = maxDouble2;
117 else
118 dfj = -maxDouble2;
119 }
120
122 // if ( ! ( dfj > 0) && ! ( dfj <= 0 ) )
123 // std::cout << " nan : dfj = " << dfj << " fval = " << fval << " invF = " << invFval << " grad = " << mfg[j] << " par[j] = " << par[j] << std::endl;
124
125 //std::cout << " x = " << currentPosition[0] << " par[j] = " << par[j] << " : dfj = " << dfj << " fval = " << fval << " invF = " << invFval << " grad = " << mfg[j] << " deriv = " << grad[j] << std::endl;
126
127
128 // in second derivative use Fumili approximation neglecting the term containing the
129 // second derivatives of the model function
130 for (int k = j; k < npar; ++ k) {
131 int idx = j + k*(k+1)/2;
132 if (std::fabs( mfg[k]) < minDouble ) {
133 if (mfg[k] < 0)
134 mfg[k] = -minDouble;
135 else
136 mfg[k] = minDouble;
137 }
138
139 double dfk = invFval * mfg[k];
140 // avoid that dfk*dfj are one small and one huge so I get a nan
141 // to avoid summing infinite and nan later when calculating the Hessian
142 if ( std::fabs(dfk) > maxDouble2 ) {
143 if (dfk > 0)
144 dfk = maxDouble2;
145 else
146 dfk = -maxDouble2;
147 }
148
149
150 h[idx] += dfj * dfk;
151 // if ( ( ! ( h[idx] > 0) && ! ( h[idx] <= 0 ) ) )
152 // std::cout << " nan : dfj = " << dfj << " fval = " << fval << " invF = " << invFval << " gradj = " << mfg[j]
153 // << " dfk = " << dfk << " gradk = "<< mfg[k] << " hess_jk = " << h[idx] << " par[k] = " << par[k] << std::endl;
154 }
155
156 } // end param loop
157
158 } // end points loop
159
160 // std::cout <<"\nEVALUATED GRADIENT and HESSIAN " << std::endl;
161 // for (int j = 0; j < npar; ++j) {
162 // std::cout << " j = " << j << " grad = " << grad[j] << std::endl;
163 // for (int k = j; k < npar; ++k) {
164 // std::cout << " k = " << k << " hess = " << Hessian(j,k) << " " << h[ j + k*(k+1)/2] << std::endl;
165 // }
166 // }
167
168 // set Value in base class
169 SetFCNValue( logLikelihood);
170
171}
172
173 } // namespace Minuit2
174
175} // namespace ROOT
#define h(i)
Definition: RSha256.hxx:106
double sqrt(double)
double log(double)
std::vector< double > & Hessian()
virtual const std::vector< double > & Gradient() const
Return cached Value of function Gradient estimated previously using the FumiliFCNBase::EvaluateAll me...
void SetFCNValue(double value)
const ParametricFunction * ModelFunction() const
Returns the model function used for the data.
virtual int GetNumberOfMeasurements() const
Accessor to the number of measurements used for calculating the maximum likelihood.
std::vector< double > Elements(const std::vector< double > &par) const
Evaluates the model function for the different measurement points and the Parameter values supplied.
virtual void EvaluateAll(const std::vector< double > &par)
Evaluate function Value, Gradient and Hessian using Fumili approximation, for values of parameters p ...
virtual const std::vector< double > & GetMeasurement(int Index) const
Accessor to the position of the measurement (x coordinate).
Function which has parameters.
virtual void SetParameters(const std::vector< double > &params) const
Sets the parameters of the ParametricFunction.
virtual std::vector< double > GetGradient(const std::vector< double > &x) const
Member function returning the Gradient of the function with respect to its variables (but without inc...
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21