ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
PaulTest4.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 #include "Minuit2/FCNBase.h"
12 #include "Minuit2/MnMigrad.h"
13 #include "Minuit2/MnMinos.h"
15 #include "Minuit2/MnPrint.h"
17 
18 #include <iostream>
19 #include <fstream>
20 
21 #ifdef USE_SEALBASE
22 #include "SealBase/Filename.h"
23 #include "SealBase/ShellEnvironment.h"
24 #endif
25 
26 using namespace ROOT::Minuit2;
27 
28 class PowerLawFunc {
29 
30 public:
31 
32  PowerLawFunc(double p0, double p1) : fP0(p0), fP1(p1) {}
33 
34  ~PowerLawFunc() {}
35 
36  double operator()(double x) const {
37  return p1()*exp(log(x)*p0());
38  }
39 
40  double p0() const {return fP0;}
41  double p1() const {return fP1;}
42 
43 private:
44  double fP0;
45  double fP1;
46 };
47 
48 class PowerLawChi2FCN : public FCNBase {
49 
50 public:
51 
52  PowerLawChi2FCN(const std::vector<double>& meas,
53  const std::vector<double>& pos,
54  const std::vector<double>& mvar) : fMeasurements(meas),
55  fPositions(pos),
56  fMVariances(mvar) {}
57 
58  ~PowerLawChi2FCN() {}
59 
60  double operator()(const std::vector<double>& par) const {
61  assert(par.size() == 2);
62  PowerLawFunc pl(par[0], par[1]);
63  double chi2 = 0.;
64 
65  for(unsigned int n = 0; n < fMeasurements.size(); n++) {
66  chi2 += ((pl(fPositions[n]) - fMeasurements[n])*(pl(fPositions[n]) - fMeasurements[n])/fMVariances[n]);
67  }
68 
69  return chi2;
70  }
71 
72  double Up() const {return 1.;}
73 
74 private:
75  std::vector<double> fMeasurements;
76  std::vector<double> fPositions;
77  std::vector<double> fMVariances;
78 };
79 
80 class PowerLawLogLikeFCN : public FCNBase {
81 
82 public:
83 
84  PowerLawLogLikeFCN(const std::vector<double>& meas,
85  const std::vector<double>& pos) :
86  fMeasurements(meas), fPositions(pos) {}
87 
88  ~PowerLawLogLikeFCN() {}
89 
90  double operator()(const std::vector<double>& par) const {
91  assert(par.size() == 2);
92  PowerLawFunc pl(par[0], par[1]);
93  double logsum = 0.;
94 
95  for(unsigned int n = 0; n < fMeasurements.size(); n++) {
96  double k = fMeasurements[n];
97  double mu = pl(fPositions[n]);
98  logsum += (k*log(mu) - mu);
99  }
100 
101  return -logsum;
102  }
103 
104  double Up() const {return 0.5;}
105 
106 private:
107  std::vector<double> fMeasurements;
108  std::vector<double> fPositions;
109 };
110 
111 int main() {
112 
113  std::vector<double> positions;
114  std::vector<double> measurements;
115  std::vector<double> var;
116  {
117 
118 #ifdef USE_SEALBASE
119  seal::Filename inputFile (seal::Filename ("$SEAL/src/MathLibs/Minuit/tests/MnSim/paul4.txt").substitute (seal::ShellEnvironment ()));
120  std::ifstream in(inputFile.Name() );
121 #else
122  std::ifstream in("paul4.txt");
123 #endif
124  if (!in) {
125  std::cerr << "Error opening input data file" << std::endl;
126  return 1;
127  }
128 
129 
130  double x = 0., y = 0., err = 0.;
131  while(in>>x>>y>>err) {
132  // if(err < 1.e-8) continue;
133  positions.push_back(x);
134  measurements.push_back(y);
135  var.push_back(err*err);
136  }
137  std::cout<<"size= "<<var.size()<<std::endl;
138  }
139  {
140  // create Chi2 FCN function
141  std::cout<<">>> test Chi2"<<std::endl;
142  PowerLawChi2FCN fFCN(measurements, positions, var);
143 
144  MnUserParameters upar;
145  upar.Add("p0", -2.3, 0.2);
146  upar.Add("p1", 1100., 10.);
147 
148  MnMigrad migrad(fFCN, upar);
149  std::cout<<"start migrad "<<std::endl;
150  FunctionMinimum min = migrad();
151  if(!min.IsValid()) {
152  //try with higher strategy
153  std::cout<<"FM is invalid, try with strategy = 2."<<std::endl;
154  MnMigrad migrad(fFCN, upar, 2);
155  min = migrad();
156  }
157  std::cout<<"minimum: "<<min<<std::endl;
158  }
159  {
160  std::cout<<">>> test log LikeliHood"<<std::endl;
161  // create LogLikelihood FCN function
162  PowerLawLogLikeFCN fFCN(measurements, positions);
163 
164  MnUserParameters upar;
165  upar.Add("p0", -2.1, 0.2);
166  upar.Add("p1", 1000., 10.);
167 
168  MnMigrad migrad(fFCN, upar);
169  std::cout<<"start migrad "<<std::endl;
170  FunctionMinimum min = migrad();
171  if(!min.IsValid()) {
172  //try with higher strategy
173  std::cout<<"FM is invalid, try with strategy = 2."<<std::endl;
174  MnMigrad migrad(fFCN, upar, 2);
175  min = migrad();
176  }
177  std::cout<<"minimum: "<<min<<std::endl;
178  }
179  {
180  std::cout<<">>> test Simplex"<<std::endl;
181  PowerLawChi2FCN chi2(measurements, positions, var);
182  PowerLawLogLikeFCN mlh(measurements, positions);
183 
184  MnUserParameters upar;
185  std::vector<double> par; par.push_back(-2.3); par.push_back(1100.);
186  std::vector<double> err; err.push_back(1.); err.push_back(1.);
187 
188  SimplexMinimizer simplex;
189 
190  std::cout<<"start simplex"<<std::endl;
191  FunctionMinimum min = simplex.Minimize(chi2, par, err);
192  std::cout<<"minimum: "<<min<<std::endl;
193  FunctionMinimum min2 = simplex.Minimize(mlh, par, err);
194  std::cout<<"minimum: "<<min2<<std::endl;
195  }
196  return 0;
197 }
double par[1]
Definition: unuranDistr.cxx:38
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
#define assert(cond)
Definition: unittest.h:542
bool Add(const std::string &, double, double)
Add free Parameter Name, Value, Error.
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:47
virtual FunctionMinimum Minimize(const FCNBase &, const std::vector< double > &, const std::vector< double > &, unsigned int stra=1, unsigned int maxfcn=0, double toler=0.1) const
tuple pl
Definition: first.py:10
static double p1(double t, double a, double b)
int main()
Definition: PaulTest4.cxx:111
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
TRObject operator()(const T1 &t1) const
Double_t y[n]
Definition: legend1.C:17
Class implementing the required methods for a minimization using Simplex.
double exp(double)
const Int_t n
Definition: legend1.C:16
double log(double)