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 
16 namespace ROOT {
17 
18 namespace Minuit2 {
19 
20 std::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 
45 const 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 
57 void 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
ROOT::Minuit2::FumiliFCNBase::SetFCNValue
void SetFCNValue(double value)
Definition: FumiliFCNBase.h:135
ROOT::Minuit2::ParametricFunction::GetGradient
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...
Definition: ParametricFunction.cxx:24
ROOT::Minuit2::ParametricFunction::SetParameters
virtual void SetParameters(const std::vector< double > &params) const
Sets the parameters of the ParametricFunction.
Definition: ParametricFunction.h:78
ROOT::Minuit2::FumiliFCNBase::Gradient
virtual const std::vector< double > & Gradient() const
Return cached Value of function Gradient estimated previously using the FumiliFCNBase::EvaluateAll me...
Definition: FumiliFCNBase.h:98
FumiliStandardMaximumLikelihoodFCN.h
ROOT::Minuit2::ParametricFunction
Function which has parameters.
Definition: ParametricFunction.h:43
log
double log(double)
ROOT::Minuit2::FumiliStandardMaximumLikelihoodFCN::EvaluateAll
virtual void EvaluateAll(const std::vector< double > &par)
Evaluate function Value, Gradient and Hessian using Fumili approximation, for values of parameters p ...
Definition: FumiliStandardMaximumLikelihoodFCN.cxx:57
ROOT::Minuit2::FumiliStandardMaximumLikelihoodFCN::Elements
std::vector< double > Elements(const std::vector< double > &par) const
Evaluates the model function for the different measurement points and the Parameter values supplied.
Definition: FumiliStandardMaximumLikelihoodFCN.cxx:20
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
h
#define h(i)
Definition: RSha256.hxx:106
ROOT::Minuit2::FumiliFCNBase::Hessian
std::vector< double > & Hessian()
Definition: FumiliFCNBase.h:139
ROOT::Minuit2::FumiliStandardMaximumLikelihoodFCN::GetMeasurement
virtual const std::vector< double > & GetMeasurement(int Index) const
Accessor to the position of the measurement (x coordinate).
Definition: FumiliStandardMaximumLikelihoodFCN.cxx:45
sqrt
double sqrt(double)
ROOT::Minuit2::FumiliMaximumLikelihoodFCN::ModelFunction
const ParametricFunction * ModelFunction() const
Returns the model function used for the data.
Definition: FumiliMaximumLikelihoodFCN.h:71
ROOT::Minuit2::FumiliStandardMaximumLikelihoodFCN::GetNumberOfMeasurements
virtual int GetNumberOfMeasurements() const
Accessor to the number of measurements used for calculating the maximum likelihood.
Definition: FumiliStandardMaximumLikelihoodFCN.cxx:51
ROOT::Minuit2::FumiliStandardMaximumLikelihoodFCN::fPositions
std::vector< std::vector< double > > fPositions
Definition: FumiliStandardMaximumLikelihoodFCN.h:129
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4