ROOT  6.06/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:823
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:6174
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3159
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
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:56
int testUnuran(TUnuran &unr, double &time, TH1 *h1, const TH1 *href, bool weightHist=false)
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:107
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
int testBinomial()
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:45
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:659
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:123
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:382
int testPoisson()
virtual void SetSeed(UInt_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
THist< 1, double > TH1D
Definition: THist.h:314
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:6669
bool useRandomSeed
int unuranDiscrete()
int testRootPoisson(double mu, double &time, TH1 *h2)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Double_t x[n]
Definition: legend1.C:17
int n
int icanv
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
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.
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:2925
TH1F * h1
Definition: legend1.C:5
virtual Double_t Chi2Test(const TH1 *h2, Option_t *option="UU", Double_t *res=0) const
chi^{2} test for comparing weighted and unweighted histograms
Definition: TH1.cxx:1850
virtual void SetBinError(Int_t bin, Double_t error)
see convention for numbering bins in TH1::GetBin
Definition: TH1.cxx:8528
virtual void SetLineColor(Color_t lcolor)
Definition: TAttLine.h:54
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2878
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:8543
TMarker * m
Definition: textangle.C:8
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7354
Double_t E()
Definition: TMath.h:54
const std::string & MethodName() const
used Unuran method
Definition: TUnuran.h:239
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
double poisson_pmf(double *x, double *p)
The Canvas class.
Definition: TCanvas.h:48
TCanvas * c1
int main(int argc, char **argv)
TUnuran class.
Definition: TUnuran.h:81
double f(double x)
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:80
double binomial_pmf(double *x, double *p)
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
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)
Automatic pad generation by division.
Definition: TPad.cxx:1077
double f2(const double *x)
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:149
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8350
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
Definition: Rtypes.h:61
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:45
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:431
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:492
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2179
int testProbVector()
virtual void Close(Option_t *option="")
Close a file.
Definition: TFile.cxx:898
Stopwatch class.
Definition: TStopwatch.h:30