Logo 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
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
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
std::vector< double > & Hessian()
virtual const std::vector< double > & Gradient() const
Return cached Value of function Gradient estimated previously using the FumiliFCNBase::EvaluateAll me...
Definition: FumiliFCNBase.h:98
void SetFCNValue(double value)
const ParametricFunction * ModelFunction() const
Returns the model function used for the data.
std::vector< double > Elements(const std::vector< double > &par) const override
Evaluates the model function for the different measurement points and the Parameter values supplied.
int GetNumberOfMeasurements() const override
Accessor to the number of measurements used for calculating the maximum likelihood.
const std::vector< double > & GetMeasurement(int Index) const override
Accessor to the position of the measurement (x coordinate).
void EvaluateAll(const std::vector< double > &par) override
Evaluate function Value, Gradient and Hessian using Fumili approximation, for values of parameters p ...
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...
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1765
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.