Logo ROOT   6.08/07
Reference Guide
DemoGaussSim.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 "GaussFcn.h"
11 #include "GaussDataGen.h"
14 #include "Minuit2/MnPrint.h"
15 #include "Minuit2/MnMigrad.h"
16 #include "Minuit2/MnMinos.h"
17 #include "Minuit2/MnContours.h"
18 #include "Minuit2/MnPlot.h"
19 #include "Minuit2/MinosError.h"
20 #include "Minuit2/ContoursError.h"
21 
22 #include <iostream>
23 
24 using namespace ROOT::Minuit2;
25 
26 int main() {
27 
28  // generate the data (100 data points)
29  GaussDataGen gdg(100);
30 
31  std::vector<double> pos = gdg.Positions();
32  std::vector<double> meas = gdg.Measurements();
33  std::vector<double> var = gdg.Variances();
34 
35  // create FCN function
36  GaussFcn fFCN(meas, pos, var);
37 
38  // create initial starting values for parameters
39  double x = 0.;
40  double x2 = 0.;
41  double norm = 0.;
42  double dx = pos[1]-pos[0];
43  double area = 0.;
44  for(unsigned int i = 0; i < meas.size(); i++) {
45  norm += meas[i];
46  x += (meas[i]*pos[i]);
47  x2 += (meas[i]*pos[i]*pos[i]);
48  area += dx*meas[i];
49  }
50  double mean = x/norm;
51  double rms2 = x2/norm - mean*mean;
52  double rms = rms2 > 0. ? sqrt(rms2) : 1.;
53 
54  {
55  // demonstrate minimal required interface for minimization
56  // create Minuit parameters without names
57 
58  // starting values for parameters
59  std::vector<double> init_par;
60  init_par.push_back(mean);
61  init_par.push_back(rms);
62  init_par.push_back(area);
63 
64  // starting values for initial uncertainties
65  std::vector<double> init_err;
66  init_err.push_back(0.1);
67  init_err.push_back(0.1);
68  init_err.push_back(0.1);
69 
70  // create minimizer (default constructor)
71  VariableMetricMinimizer fMinimizer;
72 
73  // Minimize
74  FunctionMinimum min = fMinimizer.Minimize(fFCN, init_par, init_err);
75 
76  // output
77  std::cout<<"minimum: "<<min<<std::endl;
78  }
79 
80  {
81  // demonstrate standard minimization using MIGRAD
82  // create Minuit parameters with names
83  MnUserParameters upar;
84  upar.Add("mean", mean, 0.1);
85  upar.Add("sigma", rms, 0.1);
86  upar.Add("area", area, 0.1);
87 
88  // create MIGRAD minimizer
89  MnMigrad migrad(fFCN, upar);
90 
91  // Minimize
92  FunctionMinimum min = migrad();
93 
94  // output
95  std::cout<<"minimum: "<<min<<std::endl;
96  }
97 
98  {
99  // demonstrate full interaction with parameters over subsequent
100  // minimizations
101 
102  // create Minuit parameters with names
103  MnUserParameters upar;
104  upar.Add("mean", mean, 0.1);
105  upar.Add("sigma", rms, 0.1);
106  upar.Add("area", area, 0.1);
107 
108  // access Parameter by Name to set limits...
109  upar.SetLimits("mean", mean-0.01, mean+0.01);
110 
111  // ... or access Parameter by Index
112  upar.SetLimits(1, rms-0.1, rms+0.1);
113 
114  // create Migrad minimizer
115  MnMigrad migrad(fFCN, upar);
116 
117  // Fix a Parameter...
118  migrad.Fix("mean");
119 
120  // ... and Minimize
121  FunctionMinimum min = migrad();
122 
123  // output
124  std::cout<<"minimum: "<<min<<std::endl;
125 
126  // Release a Parameter...
127  migrad.Release("mean");
128 
129  // ... and Fix another one
130  migrad.Fix(1);
131 
132  // and Minimize again
133  FunctionMinimum min1 = migrad();
134 
135  // output
136  std::cout<<"minimum1: "<<min1<<std::endl;
137 
138  // Release the Parameter...
139  migrad.Release(1);
140 
141  // ... and Minimize with all three parameters (still with limits!)
142  FunctionMinimum min2 = migrad();
143 
144  // output
145  std::cout<<"minimum2: "<<min2<<std::endl;
146 
147  // remove all limits on parameters...
148  migrad.RemoveLimits("mean");
149  migrad.RemoveLimits("sigma");
150 
151  // ... and Minimize again with all three parameters (now without limits!)
152  FunctionMinimum min3 = migrad();
153 
154  // output
155  std::cout<<"minimum3: "<<min3<<std::endl;
156  }
157 
158  {
159  // test single sided limits
160  MnUserParameters upar;
161  upar.Add("mean", mean, 0.1);
162  upar.Add("sigma", rms-1., 0.1);
163  upar.Add("area", area, 0.1);
164 
165  // test Lower limits
166  upar.SetLowerLimit("mean", mean-0.01);
167 
168  // test Upper limits
169  upar.SetUpperLimit("sigma", rms-0.5);
170 
171  // create MIGRAD minimizer
172  MnMigrad migrad(fFCN, upar);
173 
174  // ... and Minimize
175  FunctionMinimum min = migrad();
176  std::cout<<"test Lower limit minimim= "<<min<<std::endl;
177  }
178 
179  {
180  // demonstrate MINOS Error analysis
181 
182  // create Minuit parameters with names
183  MnUserParameters upar;
184  upar.Add("mean", mean, 0.1);
185  upar.Add("sigma", rms, 0.1);
186  upar.Add("area", area, 0.1);
187 
188  // create Migrad minimizer
189  MnMigrad migrad(fFCN, upar);
190 
191  // Minimize
192  FunctionMinimum min = migrad();
193 
194  // create MINOS Error factory
195  MnMinos Minos(fFCN, min);
196 
197  {
198  // 1-sigma MINOS errors (minimal interface)
199  std::pair<double,double> e0 = Minos(0);
200  std::pair<double,double> e1 = Minos(1);
201  std::pair<double,double> e2 = Minos(2);
202 
203  // output
204  std::cout<<"1-sigma Minos errors: "<<std::endl;
205  std::cout<<"par0: "<<min.UserState().Value("mean")<<" "<<e0.first<<" "<<e0.second<<std::endl;
206  std::cout<<"par1: "<<min.UserState().Value(1)<<" "<<e1.first<<" "<<e1.second<<std::endl;
207  std::cout<<"par2: "<<min.UserState().Value("area")<<" "<<e2.first<<" "<<e2.second<<std::endl;
208  }
209 
210  {
211  // 2-sigma MINOS errors (rich interface)
212  fFCN.SetErrorDef(4.);
213  MinosError e0 = Minos.Minos(0);
214  MinosError e1 = Minos.Minos(1);
215  MinosError e2 = Minos.Minos(2);
216 
217  // output
218  std::cout<<"2-sigma Minos errors: "<<std::endl;
219  std::cout<<e0<<std::endl;
220  std::cout<<e1<<std::endl;
221  std::cout<<e2<<std::endl;
222  }
223  }
224 
225  {
226  // demostrate MINOS Error analysis with limits
227 
228  // create Minuit parameters with names
229  MnUserParameters upar;
230  upar.Add("mean", mean, 0.1);
231  upar.Add("sigma", rms, 0.1);
232  upar.Add("area", area, 0.1);
233 
234  double meanLow = -50.03;
235  double rmsUp = 1.55;
236  std::cout << "sigma Limit: " << rmsUp << "\tmean limit: " << meanLow << std::endl;
237  // test Lower limits
238  upar.SetLowerLimit("mean", meanLow);
239  // test Upper limits
240  upar.SetUpperLimit("sigma", rmsUp);
241 
242  // create Migrad minimizer
243  MnMigrad migrad(fFCN, upar);
244 
245  // Minimize
246  FunctionMinimum min = migrad();
247 
248  // create MINOS Error factory
249  MnMinos Minos(fFCN, min);
250 
251  {
252  // 3-sigma MINOS errors (minimal interface)
253  fFCN.SetErrorDef(9.);
254  std::pair<double,double> e0 = Minos(0);
255  std::pair<double,double> e1 = Minos(1);
256  std::pair<double,double> e2 = Minos(2);
257 
258 
259  // output
260  std::cout<<"3-sigma Minos errors with limits: "<<std::endl;
261  std::cout.precision(16);
262  std::cout<<"par0: "<<min.UserState().Value("mean")<<" "<<e0.first<<" "<<e0.second<<std::endl;
263  std::cout<<"par1: "<<min.UserState().Value(1)<<" "<<e1.first<<" "<<e1.second<<std::endl;
264  std::cout<<"par2: "<<min.UserState().Value("area")<<" "<<e2.first<<" "<<e2.second<<std::endl;
265 
266 
267  }
268 
269  }
270 
271  {
272  // demonstrate how to use the CONTOURs
273 
274  // create Minuit parameters with names
275  MnUserParameters upar;
276  upar.Add("mean", mean, 0.1);
277  upar.Add("sigma", rms, 0.1);
278  upar.Add("area", area, 0.1);
279 
280  // create Migrad minimizer
281  MnMigrad migrad(fFCN, upar);
282 
283  // Minimize
284  FunctionMinimum min = migrad();
285 
286  // create contours factory with FCN and Minimum
287  MnContours contours(fFCN, min);
288 
289  //70% confidence level for 2 parameters Contour around the Minimum
290  // (minimal interface)
291  fFCN.SetErrorDef(2.41);
292  std::vector<std::pair<double,double> > cont = contours(0, 1, 20);
293 
294  //95% confidence level for 2 parameters Contour
295  // (rich interface)
296  fFCN.SetErrorDef(5.99);
297  ContoursError cont4 = contours.Contour(0, 1, 20);
298 
299  // plot the contours
300  MnPlot plot;
301  cont.insert(cont.end(), cont4().begin(), cont4().end());
302  plot(min.UserState().Value("mean"), min.UserState().Value("sigma"), cont);
303 
304  // print out one Contour
305  std::cout<<cont4<<std::endl;
306  }
307 
308  return 0;
309 }
int main()
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
void SetLowerLimit(unsigned int, double)
MinosError Minos(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
ask for MinosError (Lower + Upper) can be printed via std::cout
Definition: MnMinos.cxx:88
RooCmdArg Minos(Bool_t flag=kTRUE)
API class for Contours Error analysis (2-dim errors); minimization has to be done before and Minimum ...
Definition: MnContours.h:37
MnPlot produces a text-screen graphical output of (x,y) points, e.g.
Definition: MnPlot.h:26
bool Add(const std::string &, double, double)
Add free Parameter Name, Value, Error.
double sqrt(double)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
void SetLimits(unsigned int, double, double)
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
API class for Minos Error analysis (asymmetric errors); minimization has to be done before and Minimu...
Definition: MnMinos.h:34
std::vector< double > Variances() const
Definition: GaussDataGen.h:30
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
Instantiates the SeedGenerator and MinimumBuilder for Variable Metric Minimization method...
std::vector< double > Measurements() const
Definition: GaussDataGen.h:29
ContoursError Contour(unsigned int, unsigned int, unsigned int npoints=20) const
ask for one Contour ContoursError (MinosErrors + points) from number of points (>=4) and parameter in...
Definition: MnContours.cxx:39
void RemoveLimits(unsigned int)
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
std::vector< double > Positions() const
Definition: GaussDataGen.h:28
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
void SetErrorDef(double def)
add interface to set dynamically a new error definition Re-implement this function if needed...
Definition: GaussFcn.h:42
void SetUpperLimit(unsigned int, double)
const MnUserParameterState & UserState() const
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40