Logo ROOT   6.10/09
Reference Guide
unuranDistr.cxx
Go to the documentation of this file.
1 // test using 1D Distribution object interface
2 // and compare results and CPU performances using TF1::GetRandom
3 //
4 // run within ROOT (.x unuranDistr.cxx+) or pass any extra parameter in the command line to get
5 // a graphics output
6 
7 
8 #include "TStopwatch.h"
9 
10 #include "TUnuran.h"
11 #include "TUnuranContDist.h"
12 
13 #include "TH1.h"
14 #include "TF1.h"
15 
16 #include "TRandom3.h"
17 #include "TSystem.h"
18 #include "TStyle.h"
19 
20 #include "TApplication.h"
21 #include "TCanvas.h"
22 
23 #include "Math/DistFunc.h"
24 #include <cmath>
25 #include <cassert>
26 
27 #include "TError.h"
28 
29 #include <iostream>
30 
31 using std::cout;
32 using std::endl;
33 
34 int sampleSize = 5000000;
35 
36 bool useRandomSeed = false; // to use a random seed different every time
37 
38 double par[1] = {1}; // function parameters
39 
40 double norm(double *x, double *p) {
41  return ROOT::Math::normal_pdf(x[0],p[0]);
42 }
43 
44 double cdf(double *x, double *p) {
45  return ROOT::Math::normal_cdf(x[0],p[0]);
46 }
47 double cdf_trunc(double *x, double *p) {
48  double a1 = ROOT::Math::normal_cdf(p[1],p[0]);
49  double a2 = ROOT::Math::normal_cdf(p[2],p[0]);
50 
51  return ( ROOT::Math::normal_cdf(x[0],p[0]) - a1 )/(a2-a1);
52 }
53 
54 class DistTest {
55 public:
56  DistTest(TF1 * f) :
57  fCdf(f),
58  fHref(0)
59  {
60  // generate reference histo for distribution using cdf
61 
62  // make ref histo (uniform histo between 0,1
63  fHref = new TH1D("Href","uniform ref histo",100,0,1);
64  for (int i = 0; i < sampleSize; ++i)
65  fHref->Fill(gRandom->Rndm() );
66  }
67 
68  void SetCdf(TF1 *f) {
69  fCdf = f;
70  }
71 
72  int testUnuran(TUnuran & unr) {
73 
74  assert(fHref != 0);
75  assert(fCdf != 0);
76 
77  // test first the time
78  TStopwatch w;
79 
80  w.Start();
81  for (int i = 0; i < sampleSize; ++i)
82  unr.Sample();
83 
84  w.Stop();
85  double time = w.CpuTime()*1.E9/sampleSize;
86 
87 
88  TH1D htmp("htmp","gaussian generated cdf",100,0,1.);
89  // test quality (use cdf to avoid zero bins)
90  int n2 = sampleSize/100;
91  double x;
92  for (int i = 0; i<n2; ++i) {
93  x = unr.Sample();
94  htmp.Fill( fCdf->Eval(x) );
95  }
96  double prob = fHref->Chi2Test(&htmp,"UU");
97  std::cout << "Time using Unuran " << unr.MethodName() << " \t=\t " << time << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
98  if (prob < 1E-06) {
99  std::cout << "Chi2 Test failed ! " << std::endl;
100  fHref->Chi2Test(&htmp,"UUP"); // print all chi2 test info
101  return 1;
102  }
103  return 0;
104  }
105 
106  int testGetRandom(TF1 * f) {
107 
108  assert(fHref != 0);
109  assert(fCdf != 0);
110 
111  // test first the time
112 
113  TStopwatch w;
114  w.Start();
115  for (int i = 0; i < sampleSize; ++i) {
116  f->GetRandom();
117  }
118 
119  w.Stop();
120  double time = w.CpuTime()*1.E9/sampleSize;
121 
122 
123  TH1D htmp("htmp","gaussian generated cdf",100,0,1.);
124  // test quality (use cdf to avoid zero bins)
125  int n2 = sampleSize/100;
126  for (int i = 0; i<n2; ++i) {
127  double x = f->GetRandom();
128  htmp.Fill( fCdf->Eval(x) );
129  }
130  double prob = fHref->Chi2Test(&htmp,"UU");
131  std::cout << "Time using TF1::GetRandom() \t=\t " << time << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
132  if (prob < 1E-06) {
133  std::cout << "Chi2 Test failed ! " << std::endl;
134  fHref->Chi2Test(&htmp,"UUP"); // print all chi2 test info
135  return 2;
136  }
137  return 0;
138  }
139 
140 private:
141 
142  TF1 * fCdf; // cumulative distribution
143  TH1 * fHref; // reference histo
144 
145 };
146 
147 int unuranDistr() {
148 
149  int iret = 0;
150 
151  // switch off printing of info messages from chi2 test
152  gErrorIgnoreLevel = 1001;
153 
154  // check if using a random seed
155  if (useRandomSeed) gRandom->SetSeed(0);
156 
157 
158  gSystem->Load("libUnuran");
159 
160  // simple test of unuran
161 
162 
163  TH1D * h1 = new TH1D("h1","gaussian distribution",100,-10,10);
164  TH1D * h2 = new TH1D("h2","gaussian distribution",100,-10,10);
165 
166  TH1D * h1u = new TH1D("h1u","test gaussian dist",100,0,1);
167  TH1D * h2u = new TH1D("h2u","test gaussian dist",100,0,1);
168 
169 
170 
171  TF1 * f = new TF1("n",norm,-10,10,1);
172  f->SetParameters(par);
173 
174  TF1 * fc = new TF1("c",cdf,0,1,1);
175  fc->SetParameters(par);
176 
177 
178  // tester class
179  DistTest t(fc);
180 
181 
182  // create unuran 1D distribution
184 
185  // create unuran class
186  TUnuran unr(gRandom,2);
187 
188  // test all unuran methods
189  bool ret = false;
190 
191 
192  // arou: needs pdf and dpdf (estimated numerically in dist class)
193  ret = unr.Init(dist,"arou");
194  if (!ret) {
195  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
196  iret = -1;
197  }
198  else
199  iret |= t.testUnuran(unr);
200 
201  // nrou (needs only pdf , mode is an option)
202  ret = unr.Init(dist,"nrou");
203  if (!ret) {
204  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
205  iret = -2;
206  }
207  else
208  iret |= t.testUnuran(unr);
209 
210 
211  // tdr: needs pdf and dpdf (estimated numerically in dist class)
212  ret = unr.Init(dist,"tdr");
213  if (!ret) {
214  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
215  iret = -3;
216  }
217  else
218  iret |= t.testUnuran(unr);
219 
220 
221  dist.SetCdf(fc);
222 
223  // hinv (needs cdf , pdf and dpdf are optionally)
224  ret = unr.Init(dist,"hinv");
225  if (!ret) {
226  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
227  iret = -4;
228  }
229  else
230  iret |= t.testUnuran(unr);
231 
232 
233  // ninv (needs cdf, pdf is an option)
234  ret = unr.Init(dist,"ninv");
235  sampleSize /= 10; // method is too slow
236  if (!ret) {
237  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
238  iret = -5;
239  }
240  else
241  iret |= t.testUnuran(unr);
242 
243 
244  dist.SetMode( 0.0 );
245  dist.SetPdfArea(1. );
246 
247  // srou (need pdf mode sand area)
248  ret = unr.Init(dist,"srou");
249  if (!ret) {
250  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
251  iret = -6;
252  }
253  else
254  iret |= t.testUnuran(unr);
255 
256  // srou (need pdf mode sand area)
257  ret = unr.Init(dist,"ssr");
258  if (!ret) {
259  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
260  iret = -7;
261  }
262  else
263  iret |= t.testUnuran(unr);
264 
265  sampleSize *= 10;
266  ret = unr.Init(dist,"utdr");
267  if (!ret) {
268  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
269  iret = -8;
270  }
271  else
272  iret |= t.testUnuran(unr);
273 
274 
275  // test with TF1::GetRandom
276  std::cout << "\n";
277  iret |= t.testGetRandom(f);
278 
279  std::cout << "\nTest truncated distribution:\n\n";
280 
281  dist.SetDomain(-1.,1.);
282  // change cdf for tester
283  TF1 * fc2 = new TF1("fc2",cdf_trunc,-1,1,3);
284  fc2->SetParameters(par[0],-1,1.);
285  t.SetCdf(fc2);
286 
287  ret = unr.Init(dist,"auto");
288  if (!ret) {
289  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
290  iret = -10;
291  }
292  else
293  iret |= t.testUnuran(unr);
294 
295  f->SetRange(-1,1);
296  iret |= t.testGetRandom(f);
297 
298 
299 
300  TCanvas * c1 = new TCanvas("c1_unuran1D","Onedimensional distribution",10,10,1000,1000);
301  c1->Divide(2,2);
302  gStyle->SetOptFit();
303 
304  // remove the domain
305  dist.SetDomain(0,0);
306  //f->SetRange(1,-1);
307  f->SetRange(-10,10); // set verly low and large values is enough
308  // show now some plots
309  ret = unr.Init(dist,"auto");
310  if (!ret) {
311  std::cerr << "Error initializing unuran with method " << unr.MethodName() << std::endl;
312  iret = -20;
313  }
314  int n2 = sampleSize/10;
315  for (int i = 0; i < n2; ++i) {
316  double x1 = unr.Sample();
317  h1->Fill( x1 );
318  h1u->Fill( fc->Eval( x1 ) );
319  double x2 = f->GetRandom();
320  h2->Fill( x2 );
321  h2u->Fill( fc->Eval( x2 ) );
322  }
323 
324  c1->cd(1);
325  h1->Draw();
326  h1->Fit("gaus","Q");
327  std::cout << "\nFit result on data with unuran: " << std::endl;
328  TF1 * f1 = h1->GetFunction("gaus");
329  std::cout << "Gaus Fit chi2 = " << f1->GetChisquare() << " ndf = " << f1->GetNDF() << std::endl;
330  std::cout << "Gaus Fit Prob = " << f1->GetProb() << std::endl;
331  if ( f1->GetProb() < 1.E-6) iret = 11;
332  c1->cd(2);
333  h1u->Draw("Error");
334  h1u->Fit("pol0","Q");
335  TF1 * f1u = h1u->GetFunction("pol0");
336  std::cout << "CDF Fit chi2 = " << f1u->GetChisquare() << " ndf = " << f1u->GetNDF() << std::endl;
337  std::cout << "CDF Fit Prob = " << f1u->GetProb() << std::endl;
338  if ( f1u->GetProb() < 1.E-6) iret = 12;
339 
340 
341 
342 
343  c1->cd(3);
344  h2->Draw("same");
345  h2->Fit("gaus","Q");
346  std::cout << "\nFit result on data with GetRandom: " << std::endl;
347  f1 = h2->GetFunction("gaus");
348  std::cout << "Gaus Fit chi2 = " << f1->GetChisquare() << " ndf = " << f1->GetNDF() << std::endl;
349  std::cout << "Gaus Fit Prob = " << f1->GetProb() << std::endl;
350  if ( f1->GetProb() < 1.E-6) iret = 13;
351  c1->cd(4);
352 
353  h2u->Draw("Error");
354  h2u->Fit("pol0","Q");
355  f1 = h2u->GetFunction("pol0");
356  std::cout << "Fit chi2 = " << f1->GetChisquare() << " ndf = " << f1->GetNDF() << std::endl;
357  std::cout << "Fit Prob = " << f1->GetProb() << std::endl;
358  if ( f1->GetProb() < 1.E-6) iret = 14;
359 
360 
361  std::cout << "Chi2 test Gaussian histograms :\t";
362  h1->Chi2Test(h2,"UUP");
363  std::cout << "Chi2 test Uniform histograms :\t";
364  h1u->Chi2Test(h2u,"UUP");
365 
366  if (iret != 0)
367  std::cerr <<"\n\nUnuRan Continous Distribution Test:\t Failed !!!!!!!\n" << std::endl;
368  else
369  std::cerr << "\n\nUnuRan Continous Distribution Test:\t OK\n" << std::endl;
370  return iret;
371 
372 
373  return iret;
374 
375 }
376 
377 #ifndef __CLING__
378 int main(int argc, char **argv)
379 {
380  int iret = 0;
381  if (argc > 1) {
382  TApplication theApp("App",&argc,argv);
383  iret = unuranDistr();
384  theApp.Run();
385  }
386  else
387  iret = unuranDistr();
388  return iret;
389 }
390 #endif
double par[1]
Definition: unuranDistr.cxx:38
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3126
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:588
virtual Double_t GetRandom()
Return a random number following this function shape.
Definition: TF1.cxx:1815
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
int testUnuran(TUnuran &unr, double &time, TH1 *h1, const TH1 *href, bool weightHist=false)
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:105
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:389
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition: TH1.cxx:8163
return c1
Definition: legend1.C:41
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:679
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF1.cxx:3288
int unuranDistr()
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1825
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
void SetCdf(TF1 *cdf)
set cdf distribution.
virtual Double_t GetProb() const
Return the fit probability.
Definition: TF1.cxx:1681
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
bool useRandomSeed
Definition: unuranDistr.cxx:36
virtual Int_t GetNDF() const
Return the number of degrees of freedom in the fit the fNDF parameter has been previously computed du...
Definition: TF1.cxx:1616
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:1956
int sampleSize
Definition: unuranDistr.cxx:34
static const double x2[5]
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
virtual Double_t Chi2Test(const TH1 *h2, Option_t *option="UU", Double_t *res=0) const
test for comparing weighted and unweighted histograms
Definition: TH1.cxx:1830
double cdf(double *x, double *p)
Definition: unuranDistr.cxx:44
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
double cdf_trunc(double *x, double *p)
Definition: unuranDistr.cxx:47
int testGetRandom(TF2 *f, TH1 *h1, const TH2 *href=0)
TH1F * h1
Definition: legend1.C:5
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:512
R__EXTERN TSystem * gSystem
Definition: TSystem.h:539
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1219
int main(int argc, char **argv)
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:594
constexpr Double_t E()
Definition: TMath.h:74
The Canvas class.
Definition: TCanvas.h:31
Double_t GetChisquare() const
Definition: TF1.h:398
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Evaluate this function.
Definition: TF1.cxx:1197
void SetDomain(double xmin, double xmax)
Set the distribution domain.
TUnuran class.
Definition: TUnuran.h:79
static const double x1[5]
double f(double x)
The TH1 histogram class.
Definition: TH1.h:56
void SetPdfArea(double area)
set the area below the pdf
const std::string & MethodName() const
used Unuran method
Definition: TUnuran.h:237
TUnuranContDist class describing one dimensional continuous distribution.
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1135
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:79
1-Dim function class
Definition: TF1.h:150
TF1 * f1
Definition: legend1.C:11
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:310
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
void SetMode(double mode)
set the distribution mode (x position of its maximum)
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3564
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
Stopwatch class.
Definition: TStopwatch.h:28