ROOT   Reference Guide
Loading...
Searching...
No Matches
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
11
12#include <vector>
13#include <cmath>
14#include <limits>
15
16namespace ROOT {
17
18namespace Minuit2 {
19
20std::vector<double> FumiliStandardMaximumLikelihoodFCN::Elements(const std::vector<double> &par) const
21{
22
23 // calculate likelihood element f(i) = pdf(x(i))
24 std::vector<double> result;
25 double tmp1 = 0.0;
26 unsigned int fPositionsSize = fPositions.size();
27
28 for (unsigned int i = 0; i < fPositionsSize; i++) {
29
30 const std::vector<double> &currentPosition = fPositions[i];
31
32 // The commented line is the object-oriented way to do it
33 // but it is faster to do a single function call...
34 //(*(this->getModelFunction())).SetParameters(par);
35 tmp1 = (*(this->ModelFunction()))(par, currentPosition);
36
37 // std::cout << " i = " << i << " " << currentPosition[0] << " " << tmp1 << std::endl;
38
39 result.push_back(tmp1);
40 }
41
42 return result;
43}
44
45const std::vector<double> &FumiliStandardMaximumLikelihoodFCN::GetMeasurement(int index) const
46{
47 // Return x(i).
48 return fPositions[index];
49}
50
52{
53 // return size of positions (coordinates).
54 return fPositions.size();
55}
56
57void FumiliStandardMaximumLikelihoodFCN::EvaluateAll(const std::vector<double> &par)
58{
59 // Evaluate in one loop likelihood value, gradient and hessian
60
61 const double minDouble = 8.0 * std::numeric_limits<double>::min();
62 const double minDouble2 = std::sqrt(minDouble);
63 const double maxDouble2 = 1.0 / minDouble2;
64 // loop on the measurements
65
66 int nmeas = GetNumberOfMeasurements();
67 std::vector<double> &grad = Gradient();
68 std::vector<double> &h = Hessian();
69 int npar = par.size();
70 double logLikelihood = 0;
71 grad.resize(npar);
72 h.resize(static_cast<unsigned int>(0.5 * npar * (npar + 1)));
73 grad.assign(npar, 0.0);
74 h.assign(static_cast<unsigned int>(0.5 * npar * (npar + 1)), 0.0);
75
76 const ParametricFunction &modelFunc = *ModelFunction();
77
78 for (int i = 0; i < nmeas; ++i) {
79
80 // work for one-dimensional points
81 const std::vector<double> &currentPosition = fPositions[i];
82 modelFunc.SetParameters(currentPosition);
83 double fval = modelFunc(par);
84 if (fval < minDouble)
85 fval = minDouble; // to avoid getting infinity and nan's
86 logLikelihood -= std::log(fval);
87 double invFval = 1.0 / fval;
88 // this method should return a reference
89 std::vector<double> mfg = modelFunc.GetGradient(par);
90
91 // calc derivatives
92
93 for (int j = 0; j < npar; ++j) {
94 if (std::fabs(mfg[j]) < minDouble) {
95 // std::cout << "SMALL values: grad = " << mfg[j] << " " << minDouble << " f(x) = " << fval
96 // << " params " << j << " p0 = " << par[0] << " p1 = " << par[1] << std::endl;
97 if (mfg[j] < 0)
98 mfg[j] = -minDouble;
99 else
100 mfg[j] = minDouble;
101 }
102
103 double dfj = invFval * mfg[j];
104 // to avoid summing infinite and nan later when calculating the Hessian
105 if (std::fabs(dfj) > maxDouble2) {
106 if (dfj > 0)
107 dfj = maxDouble2;
108 else
109 dfj = -maxDouble2;
110 }
111
112 grad[j] -= dfj;
113 // if ( ! ( dfj > 0) && ! ( dfj <= 0 ) )
114 // std::cout << " nan : dfj = " << dfj << " fval = " << fval << " invF = " << invFval << " grad = " << mfg[j]
115 // << " par[j] = " << par[j] << std::endl;
116
117 // std::cout << " x = " << currentPosition[0] << " par[j] = " << par[j] << " : dfj = " << dfj << " fval = "
118 // << fval << " invF = " << invFval << " grad = " << mfg[j] << " deriv = " << grad[j] << std::endl;
119
120 // in second derivative use Fumili approximation neglecting the term containing the
121 // second derivatives of the model function
122 for (int k = j; k < npar; ++k) {
123 int idx = j + k * (k + 1) / 2;
124 if (std::fabs(mfg[k]) < minDouble) {
125 if (mfg[k] < 0)
126 mfg[k] = -minDouble;
127 else
128 mfg[k] = minDouble;
129 }
130
131 double dfk = invFval * mfg[k];
132 // avoid that dfk*dfj are one small and one huge so I get a nan
133 // to avoid summing infinite and nan later when calculating the Hessian
134 if (std::fabs(dfk) > maxDouble2) {
135 if (dfk > 0)
136 dfk = maxDouble2;
137 else
138 dfk = -maxDouble2;
139 }
140
141 h[idx] += dfj * dfk;
142 // if ( ( ! ( h[idx] > 0) && ! ( h[idx] <= 0 ) ) )
143 // std::cout << " nan : dfj = " << dfj << " fval = " << fval << " invF = " << invFval << " gradj = " <<
144 // mfg[j]
145 // << " dfk = " << dfk << " gradk = "<< mfg[k] << " hess_jk = " << h[idx] << " par[k] = " << par[k] <<
146 // std::endl;
147 }
148
149 } // end param loop
150
151 } // end points loop
152
153 // std::cout <<"\nEVALUATED GRADIENT and HESSIAN " << std::endl;
154 // for (int j = 0; j < npar; ++j) {
155 // std::cout << " j = " << j << " grad = " << grad[j] << std::endl;
156 // for (int k = j; k < npar; ++k) {
157 // std::cout << " k = " << k << " hess = " << Hessian(j,k) << " " << h[ j + k*(k+1)/2] << std::endl;
158 // }
159 // }
160
161 // set Value in base class
162 SetFCNValue(logLikelihood);
163}
164
165} // namespace Minuit2
166
167} // namespace ROOT
#define h(i)
Definition RSha256.hxx:106
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...
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...