ROOT   6.10/09 Reference Guide
unuranDiscrete.cxx
Go to the documentation of this file.
1 // test using discrete distribution.
2 // Generate numbers from a given probability vector or from a discrete distribution like
3 // the Poisson distribution.
4 // Compare also the Unuran method for generating Poisson numbers with TRandom::Poisson
5 //
6 // run within ROOT (.x unuranDiscrete.cxx+) or pass any extra parameter in the command line to get
7 // a graphics output (./unuranDiscrete 1 )
8
9
10 #include "TUnuran.h"
11 #include "TUnuranDiscrDist.h"
12
13 #include "TH1.h"
14 #include "TMath.h"
15 #include "TF1.h"
16 #include "TRandom3.h"
17 #include "TApplication.h"
18 #include "TCanvas.h"
19 #include "TStopwatch.h"
20 #include "TError.h"
21
22
23 #ifdef WRITE_OUTPUT
24 #include "TFile.h"
25 #include "TGraph.h"
26 #endif
27
28 #include "Math/Functor.h"
29 #include "Math/DistFunc.h"
30 #include "Math/Util.h"
31
32 #include <iostream>
33 #include <iomanip>
34
35 //#define DEBUG
36 #ifdef DEBUG
37 int n = 100;
38 #else
39 int n = 1000000;
40 #endif
42 int icanv = 1;
43
44 bool useRandomSeed = false; // to use a random seed different every time
45
46 double poisson_pmf(double * x, double * p) {
47
48  double y = ROOT::Math::poisson_pdf(int(x[0]),p[0]);
49 // std::cout << x[0] << " f(x) = " << y << std::endl;
50  return y;
51 }
52 double binomial_pmf(double * x, double * p) {
53
54  double y = ROOT::Math::binomial_pdf(static_cast<unsigned int>(x[0]),p[1], static_cast<unsigned int>(p[0]));
55 // std::cout << x[0] << " f(x) = " << y << std::endl;
56  return y;
57 }
58
59 int testUnuran(TUnuran & unr, double & time, TH1 * h1, const TH1 * href, bool weightHist = false ) {
60
61
62  // test first the time
63  TStopwatch w;
64
65  w.Start();
66  for (int i = 0; i < n; ++i)
67  unr.SampleDiscr();
68
69  w.Stop();
70  time = w.CpuTime()*1.E9/n;
71
72
73  // test quality (use cdf to avoid zero bins)
74  h1->Reset();
75  int n2 = n/10;
76  int x;
77  for (int i = 0; i<n2; ++i) {
78  x = unr.SampleDiscr();
79  h1->Fill( double( x) );
80  }
81  double prob;
82  if (weightHist)
83  prob = h1->Chi2Test(href,"UW");
84  else
85  prob = h1->Chi2Test(href,"UU");
86  std::string s = "Time using Unuran " + unr.MethodName();
87  std::cout << std::left << std::setw(40) << s << "\t=\t " << time << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
88  if (prob < 1E-06) {
89  std::cout << "Chi2 Test failed for method " << unr.MethodName() << std::endl;
90  if (weightHist)
91  h1->Chi2Test(href,"UWP"); // print all chi2 test info
92  else
93  h1->Chi2Test(href,"UUP"); // print all chi2 test info
94
95  return 1;
96  }
97  return 0;
98 }
99
100 int testRootPoisson(double mu, double &time, TH1 * h2) {
101  TStopwatch w;
102
103  w.Start();
104  for (int i = 0; i < n; ++i)
105  gRandom->Poisson(mu);
106
107  w.Stop();
108  time = w.CpuTime()*1.E9/n;
109
110  // make ref histo
111  int n2 = n/10;
112  for (int i = 0; i< n2; ++i) {
113  h2->Fill ( gRandom->Poisson(mu) );
114  }
115
116  std::string s = "Time using TRandom::Poisson";
117  std::cout << std::left << std::setw(40) << s << "\t=\t " << time << std::endl;
118  return 0;
119 }
120
121 int testRootBinomial(int m, double p, double &time, TH1 * h2) {
122
123  TStopwatch w;
124
125  w.Start();
126  for (int i = 0; i < n; ++i)
127  gRandom->Binomial(m,p);
128
129  w.Stop();
130  time = w.CpuTime()*1.E9/n;
131
132  // make ref histo
133  int n2 = n/10;
134  for (int i = 0; i< n2; ++i) {
135  h2->Fill ( gRandom->Binomial(m,p) );
136  }
137
138  std::string s = "Time using TRandom::Binomial";
139  std::cout << std::left << std::setw(40) << s << "\t=\t " << time << std::endl;
140  return 0;
141 }
142
144
145  int iret = 0;
146
147  TH1D * h0 = new TH1D("h0","reference prob",10,-0.5,9.5);
148  TH1D * h1 = new TH1D("h1","UNURAN gen prob",10,-0.5,9.5);
149
150
151  double p[10] = {1.,2.,3.,5.,3.,2.,1.,0.5,0.3,0.5 };
152  for (int i = 0; i< 10; ++i) {
153  h0->SetBinContent(i+1,p[i]);
154  h0->SetBinError(i+1,0.);
155  }
156  double sum = h0->GetSumOfWeights();
157  std::cout << " prob sum = " << sum << std::endl;
158
159  TUnuran unr(gRandom,2);
160
161  TUnuranDiscrDist dist(p,p+10);
162
163 // TUnuranDiscrDist fDDist = dist;
164 // std::cout << fDDist.ProbVec().size() << std::endl;
165
166  std::cout << "Test generation with a PV :\n\n";
167
168  double time;
169
170  bool ret = unr.Init(dist,"method=dau");
171  if (!ret) iret = -1;
172  iret |= testUnuran(unr,time,h1,h0,true);
173
174  ret = unr.Init(dist,"method=dgt");
175  if (!ret) iret = -1;
176  iret |= testUnuran(unr,time,h1,h0,true);
177
178  // dss require the sum
179  dist.SetProbSum(sum);
180  ret = unr.Init(dist,"method=dss");
181  if (!ret) iret = -1;
182  iret |= testUnuran(unr,time,h1,h0,true);
183
184
185
186
187
188  c1->cd(icanv++);
189  h1->Sumw2();
190  h1->Scale( h0->GetSumOfWeights()/(h1->GetSumOfWeights() ) );
191  h0->Draw();
192  h1->SetLineColor(kBlue);
193  h1->Draw("Esame");
194  c1->Update();
195
196  return iret;
197 }
198
199 int testPoisson() {
200  // test with a function
201  // Poisson distribution
202
203  int iret = 0;
204  std::cout << "\nTest generation with a Probability function (Poisson) :\n\n";
205
206  TF1 * f = new TF1("f",poisson_pmf,1,0,1);
207
208  // loop on mu values for Nmu times
209
210  const int Nmu = 5;
211  double muVal[Nmu] = {5,10,20,50,100};
212  double tR[Nmu];
213  double tU[Nmu];
214  double tUdari[Nmu];
215  double tUdsrou[Nmu];
216
217  TUnuran unr(gRandom,2);
218
219  for (int imu = 0; imu < Nmu; ++imu) {
220
221  const double mu = muVal[imu];
222
223  int nmax = static_cast<int>(mu*3);
224  TH1D * h2 = new TH1D("h2","reference Poisson prob",nmax,0,nmax);
225  TH1D * h3 = new TH1D("h3","UNURAN gen prob",nmax,0,nmax);
226  TH1D * h4 = new TH1D("h4","UNURAN gen prob 2",nmax,0,nmax);
227
228
229  std::cout << "\nPoisson mu = " << mu << std::endl << std::endl;
230
231  testRootPoisson(mu,tR[imu],h2);
232
233  // test UNURAN with standard method
234  bool ret = unr.InitPoisson(mu);
235  if (!ret) iret = -1;
236  testUnuran(unr,tU[imu],h3,h2);
237
238
239  // test changing all the time the mu
240  // use re-init for a fast re-initialization
241  TStopwatch w;
242  unr.InitPoisson(mu,"dstd");
243  double p[1]; p[0] = mu;
244  w.Start();
245  for (int i = 0; i < n; ++i) {
246  unr.ReInitDiscrDist(1,p);
247  int k = unr.SampleDiscr();
248  if (n % 10 == 0) h4->Fill(k);
249  }
250  w.Stop();
251  double time = w.CpuTime()*1.E9/n;
252  double prob = h2->Chi2Test(h4,"UU");
253  std::string s = "Time using Unuran w/ re-init method=" + unr.MethodName();
254  std::cout << std::left << std::setw(40) << s << "\t=\t " << time
255  << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
256
257  if (prob < 1E-06) {
258  std::cout << "Chi2 Test failed for re-init " << std::endl;
259  iret = -2;
260  }
261
262  f->SetParameter(0,mu);
263
264 #ifdef USE_FUNCTOR
267 #else
269 #endif
270
271  // dari method (needs mode and pdf sum)
272  dist2.SetMode(int(mu) );
273  dist2.SetProbSum(1. );
274  ret = unr.Init(dist2,"method=dari");
275  if (!ret) iret = -1;
276
277  iret |= testUnuran(unr,tUdari[imu],h3,h2);
278
279  // dsrou method (needs mode and sum)
280
281  dist2.SetProbSum(1);
282  ret = unr.Init(dist2,"method=dsrou");
283  if (!ret) iret = -1;
284
285  iret |= testUnuran(unr,tUdsrou[imu],h3,h2);
286
287
288  c1->cd(icanv++);
289  h2->DrawCopy();
290  h3->SetLineColor(kBlue);
291  h3->DrawCopy("Esame");
292  c1->Update();
293
294  delete h2;
295  delete h3;
296  delete h4;
297  }
298
299 #ifdef WRITE_OUTPUT
300
301  TFile * file = new TFile("unuranPoisson.root","RECREATE");
302  // create graphs with results
303  TGraph * gR = new TGraph(Nmu,muVal,tR);
304  TGraph * gU = new TGraph(Nmu,muVal,tU);
305  TGraph * gU2 = new TGraph(Nmu,muVal,tUdari);
306  TGraph * gU3 = new TGraph(Nmu,muVal,tUdsrou);
307  gR->Write();
308  gU->Write();
309  gU2->Write();
310  gU3->Write();
311  file->Close();
312
313 #endif
314
315  return iret;
316 }
317
319  // test using binomial distribution
320  int iret = 0;
321
322  std::cout << "\nTest generation with a Probability function (Binomimal) :\n\n";
323
324  TF1 * f = new TF1("f",binomial_pmf,1,0,2);
325
326  // loop on binomual values
327
328  const int NBin = 3;
329  double pVal[NBin] = {0.5,0.1,0.01};
330  double NVal[NBin] = {20,100,1000};
331  double tR[NBin];
332  double tU[NBin];
333  double tUdari[NBin];
334  double tUdsrou[NBin];
335
336
337  TUnuran unr(gRandom,2);
338
339
340  for (int ib = 0; ib < NBin; ++ib) {
341
342
343  double par[2];
344  par[0] = NVal[ib];
345  par[1] = pVal[ib];
346
347  int nmax = static_cast<int>(par[0]*par[1]*3);
348  TH1D * h2 = new TH1D("h2","reference Binomial prob",nmax,0,nmax);
349  TH1D * h3 = new TH1D("h3","UNURAN gen prob",nmax,0,nmax);
350  TH1D * h4 = new TH1D("h4","UNURAN gen prob 2",nmax,0,nmax);
351
352
353  std::cout << "\nBinomial n = " << par[0] << " " << par[1] << std::endl;
354
355  testRootBinomial(static_cast<int>(par[0]),par[1],tR[ib],h2);
356
357
358  // test UNURAN with standard method
359  bool ret = unr.InitBinomial(static_cast<int>(par[0]),par[1]);
360  if (!ret) iret = -1;
361  testUnuran(unr,tU[ib],h3,h2);
362
363
364  // test changing all the time the mu
365  // use re-init for a fast re-initialization
366
367  TStopwatch w;
368  unr.InitBinomial(static_cast<int>(par[0]), par[1],"dstd");
369  w.Start();
370  for (int i = 0; i < n; ++i) {
371  unr.ReInitDiscrDist(2,par);
372  int k = unr.SampleDiscr();
373  if (n % 10 == 0) h4->Fill(k);
374  }
375  w.Stop();
376  double time = w.CpuTime()*1.E9/n;
377  double prob = h2->Chi2Test(h4,"UU");
378  std::string s = "Time using Unuran w/ re-init method=" + unr.MethodName();
379  std::cout << std::left << std::setw(40) << s << "\t=\t " << time
380  << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
381
382  if (prob < 1E-06) {
383  std::cout << "Chi2 Test failed for re-init " << std::endl;
384  iret = -2;
385  }
386
387  // test the universal methods
388
389  f->SetParameters(par);
391
392  // dari method (needs mode and pdf sum)
393  dist.SetMode(int(par[0]*par[1]) );
394  dist.SetProbSum(1. );
395  ret = unr.Init(dist,"method=dari");
396  if (!ret) iret = -1;
397
398  iret |= testUnuran(unr,tUdari[ib],h3,h2);
399
400  // dsrou method (needs mode and sum)
401
402  ret = unr.Init(dist,"method=dsrou");
403  if (!ret) iret = -1;
404
405  iret |= testUnuran(unr,tUdsrou[ib],h3,h2);
406
407
408  c1->cd(icanv++);
409  h2->DrawCopy();
410  h3->SetLineColor(kBlue);
411  h3->DrawCopy("Esame");
412  c1->Update();
413
414  delete h2;
415  delete h3;
416  delete h4;
417  }
418
419
420  return iret;
421 }
422
424
425  int iret = 0;
426
427  c1 = new TCanvas("c1_unuranDiscr_PV","Discrete distribution from PV",10,10,800,800);
428  c1->Divide(3,3);
429
430  // switch off printing of info messages from chi2 test
431  gErrorIgnoreLevel = 1001;
432
433  // check if using a random seed
434  if (useRandomSeed) gRandom->SetSeed(0);
435
436  iret |= testProbVector();
437  iret |= testPoisson();
438  iret |= testBinomial();
439
440  if (iret != 0)
441  std::cerr <<"\n\nUnuRan Discrete Distribution Test:\t Failed !!!!!!!" << std::endl;
442  else
443  std::cerr << "\n\nUnuRan Discrete Distribution Test:\t OK" << std::endl;
444
445  return iret;
446
447 }
448
449
450 #ifndef __CINT__
451 int main(int argc, char **argv)
452 {
453  int iret = 0;
454  if (argc > 1) {
455  TApplication theApp("App",&argc,argv);
456  iret = unuranDiscrete();
457  theApp.Run();
458  }
459  else
460  iret = unuranDiscrete();
461
462  return iret;
463 }
464 #endif
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:778
double par[1]
Definition: unuranDistr.cxx:38
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:5937
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
static long int sum(long int i)
Definition: Factory.cxx:2162
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:588
void SetProbSum(double sum)
set the value of the sum of the probabilities in the given domain
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
bool InitBinomial(unsigned int ntot, double prob, const std::string &method="dstd")
Initialize method for the Binomial distribution Used to generate poisson numbers for a constant param...
Definition: TUnuran.cxx:438
virtual Int_t Binomial(Int_t ntot, Double_t prob)
Generates a random integer N according to the binomial law.
Definition: TRandom.cxx:172
bool ReInitDiscrDist(unsigned int npar, double *params)
Reinitialize UNURAN by changing the distribution parameters but mantaining same distribution and meth...
Definition: TUnuran.cxx:453
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7103
int testBinomial()
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:46
Definition: TCanvas.cxx:679
void SetMode(int mode)
set the mode of the distribution (location of maximum probability)
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
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:2898
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:382
int testPoisson()
static ROOT::R::TRInterface * gR
Definition: TRInterface.cxx:22
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:6419
bool useRandomSeed
int unuranDiscrete()
int testRootPoisson(double mu, double &time, TH1 *h2)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
int n
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
int icanv
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
bool InitPoisson(double mu, const std::string &method="dstd")
Initialize method for the Poisson distribution Used to generate poisson numbers for a constant parame...
Definition: TUnuran.cxx:423
TUnuranDiscrDist class for one dimensional discrete distribution.
TH1F * h1
Definition: legend1.C:5
virtual void SetBinError(Int_t bin, Double_t error)
See convention for numbering bins in TH1::GetBin.
Definition: TH1.cxx:8311
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
double binomial_pdf(unsigned int k, double p, unsigned int n)
Probability density function of the binomial distribution.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8325
TMarker * m
Definition: textangle.C:8
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
double poisson_pmf(double *x, double *p)
The Canvas class.
Definition: TCanvas.h:31
TCanvas * c1
int main(int argc, char **argv)
TUnuran class.
Definition: TUnuran.h:79
double f(double x)
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:56
double binomial_pmf(double *x, double *p)
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
const std::string & MethodName() const
used Unuran method
Definition: TUnuran.h:237
int testRootBinomial(int m, double p, double &time, TH1 *h2)
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)
double f2(const double *x)
Definition: file.py:1
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
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8132
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:310
Definition: Rtypes.h:56
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:578
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:362
Functor1D class for one-dimensional functions.
Definition: Functor.h:487
virtual void Update()