Logo ROOT   master
Reference Guide
unuranDemo.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_unuran
3 /// Example macro to show unuran capabilities
4 /// The results are compared with what is obtained using TRandom or TF1::GetRandom
5 /// The macro is divided in 3 parts:
6 /// - testStringAPI : show how to use string API of UNURAN to generate Gaussian random numbers
7 /// - testDistr1D : show how to pass a 1D distribution object to UNURAN to generate numbers
8 /// according to the given distribution object
9 /// - testDistrMultiDIm : show how to pass a multidimensional distribution object to UNURAN
10 ///
11 ///
12 /// To execute the macro type in:
13 ///
14 /// ~~~{.cpp}
15 /// root[0]: gSystem->Load("libMathCore");
16 /// root[0]: gSystem->Load("libUnuran");
17 /// root[0]: .x unuranDemo.C+
18 /// ~~~
19 ///
20 /// \macro_code
21 ///
22 /// \author Lorenzo Moneta
23 
24 
25 #include "TStopwatch.h"
26 
27 #include "TUnuran.h"
28 #include "TUnuranContDist.h"
29 #include "TUnuranMultiContDist.h"
30 #include "TUnuranDiscrDist.h"
31 #include "TUnuranEmpDist.h"
32 
33 #include "TH1.h"
34 #include "TH3.h"
35 #include "TF3.h"
36 #include "TMath.h"
37 
38 #include "TRandom2.h"
39 #include "TSystem.h"
40 #include "TStyle.h"
41 
42 #include "TApplication.h"
43 #include "TCanvas.h"
44 
45 #include "Math/ProbFunc.h"
46 #include "Math/DistFunc.h"
47 
48 #include <iostream>
49 #include <cassert>
50 
51 using std::cout;
52 using std::endl;
53 
54 // number of distribution generated points
55 #define NGEN 1000000
56 
57 int izone = 0;
58 TCanvas * c1 = 0;
59 
60 // test using UNURAN string interface
61 void testStringAPI() {
62 
63  TH1D * h1 = new TH1D("h1G","gaussian distribution from Unuran",100,-10,10);
64  TH1D * h2 = new TH1D("h2G","gaussian distribution from TRandom",100,-10,10);
65 
66  cout << "\nTest using UNURAN string API \n\n";
67 
68 
69  TUnuran unr;
70  if (! unr.Init( "normal()", "method=arou") ) {
71  cout << "Error initializing unuran" << endl;
72  return;
73  }
74 
75  int n = NGEN;
76  TStopwatch w;
77  w.Start();
78 
79  for (int i = 0; i < n; ++i) {
80  double x = unr.Sample();
81  h1->Fill( x );
82  }
83 
84  w.Stop();
85  cout << "Time using Unuran method " << unr.MethodName() << "\t=\t " << w.CpuTime() << endl;
86 
87 
88  // use TRandom::Gaus
89  w.Start();
90  for (int i = 0; i < n; ++i) {
91  double x = gRandom->Gaus(0,1);
92  h2->Fill( x );
93  }
94 
95  w.Stop();
96  cout << "Time using TRandom::Gaus \t=\t " << w.CpuTime() << endl;
97 
98  assert(c1 != 0);
99  c1->cd(++izone);
100  h1->Draw();
101  c1->cd(++izone);
102  h2->Draw();
103 
104 }
105 
106 
107 
108 double distr(double *x, double *p) {
109  return ROOT::Math::breitwigner_pdf(x[0],p[0],p[1]);
110 }
111 
112 double cdf(double *x, double *p) {
113  return ROOT::Math::breitwigner_cdf(x[0],p[0],p[1]);
114 }
115 
116 // test of unuran passing as input a distribution object( a BreitWigner) distribution
117 void testDistr1D() {
118 
119  cout << "\nTest 1D Continous distributions\n\n";
120 
121  TH1D * h1 = new TH1D("h1BW","Breit-Wigner distribution from Unuran",100,-10,10);
122  TH1D * h2 = new TH1D("h2BW","Breit-Wigner distribution from GetRandom",100,-10,10);
123 
124 
125 
126  TF1 * f = new TF1("distrFunc",distr,-10,10,2);
127  double par[2] = {1,0}; // values are gamma and mean
128  f->SetParameters(par);
129 
130  TF1 * fc = new TF1("cdfFunc",cdf,-10,10,2);
131  fc->SetParameters(par);
132 
133  // create Unuran 1D distribution object
135  // to use a different random number engine do:
136  TRandom2 * random = new TRandom2();
137  int logLevel = 2;
138  TUnuran unr(random,logLevel);
139 
140  // select unuran method for generating the random numbers
141  std::string method = "tdr";
142  //std::string method = "method=auto";
143  // "method=hinv"
144  // set the cdf for some methods like hinv that requires it
145  // dist.SetCdf(fc);
146 
147  //cout << "unuran method is " << method << endl;
148 
149  if (!unr.Init(dist,method) ) {
150  cout << "Error initializing unuran" << endl;
151  return;
152  }
153 
154 
155 
156  TStopwatch w;
157  w.Start();
158 
159  int n = NGEN;
160  for (int i = 0; i < n; ++i) {
161  double x = unr.Sample();
162  h1->Fill( x );
163  }
164 
165  w.Stop();
166  cout << "Time using Unuran method " << unr.MethodName() << "\t=\t " << w.CpuTime() << endl;
167 
168  w.Start();
169  for (int i = 0; i < n; ++i) {
170  double x = f->GetRandom();
171  h2->Fill( x );
172  }
173 
174  w.Stop();
175  cout << "Time using TF1::GetRandom() \t=\t " << w.CpuTime() << endl;
176 
177  c1->cd(++izone);
178  h1->Draw();
179 
180  c1->cd(++izone);
181  h2->Draw();
182 
183  std::cout << " chi2 test of UNURAN vs GetRandom generated histograms: " << std::endl;
184  h1->Chi2Test(h2,"UUP");
185 
186 }
187 
188 // 3D gaussian distribution
189 double gaus3d(double *x, double *p) {
190 
191  double sigma_x = p[0];
192  double sigma_y = p[1];
193  double sigma_z = p[2];
194  double rho = p[2];
195  double u = x[0] / sigma_x ;
196  double v = x[1] / sigma_y ;
197  double w = x[2] / sigma_z ;
198  double c = 1 - rho*rho ;
199  double result = (1 / (2 * TMath::Pi() * sigma_x * sigma_y * sigma_z * sqrt(c)))
200  * exp (-(u * u - 2 * rho * u * v + v * v + w*w) / (2 * c));
201  return result;
202 }
203 
204 // test of unuran passing as input a multi-dimension distribution object
205 void testDistrMultiDim() {
206 
207  cout << "\nTest Multidimensional distributions\n\n";
208 
209  TH3D * h1 = new TH3D("h13D","gaussian 3D distribution from Unuran",50,-10,10,50,-10,10,50,-10,10);
210  TH3D * h2 = new TH3D("h23D","gaussian 3D distribution from GetRandom",50,-10,10,50,-10,10,50,-10,10);
211 
212 
213 
214  TF3 * f = new TF3("g3d",gaus3d,-10,10,-10,10,-10,10,3);
215  double par[3] = {2,2,0.5};
216  f->SetParameters(par);
217 
218 
219 
221  TUnuran unr(gRandom);
222  //std::string method = "method=vnrou";
223 
224  //std::string method = "method=hitro;use_boundingrectangle=false ";
225  std::string method = "hitro";
226  if ( ! unr.Init(dist,method) ) {
227  cout << "Error initializing unuran" << endl;
228  return;
229  }
230 
231  TStopwatch w;
232  w.Start();
233 
234  double x[3];
235  for (int i = 0; i < NGEN; ++i) {
236  unr.SampleMulti(x);
237  h1->Fill(x[0],x[1],x[2]);
238  }
239 
240  w.Stop();
241  cout << "Time using Unuran method " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
242 
243  assert(c1 != 0);
244  c1->cd(++izone);
245  h1->Draw();
246 
247 
248  // need to set a reasonable number of points in TF1 to get acceptable quality from GetRandom to
249  int np = 200;
250  f->SetNpx(np);
251  f->SetNpy(np);
252  f->SetNpz(np);
253 
254  w.Start();
255  for (int i = 0; i < NGEN; ++i) {
256  f->GetRandom3(x[0],x[1],x[2]);
257  h2->Fill(x[0],x[1],x[2]);
258  }
259 
260  w.Stop();
261  cout << "Time using TF1::GetRandom \t\t=\t " << w.CpuTime() << endl;
262 
263 
264  c1->cd(++izone);
265  h2->Draw();
266 
267  std::cout << " chi2 test of UNURAN vs GetRandom generated histograms: " << std::endl;
268  h1->Chi2Test(h2,"UUP");
269 
270 }
271 //_____________________________________________
272 //
273 // example of discrete distributions
274 
275 double poisson(double * x, double * p) {
276  return ROOT::Math::poisson_pdf(int(x[0]),p[0]);
277 }
278 
279 void testDiscDistr() {
280 
281  cout << "\nTest Discrete distributions\n\n";
282 
283  TH1D * h1 = new TH1D("h1PS","Unuran Poisson prob",20,0,20);
284  TH1D * h2 = new TH1D("h2PS","Poisson dist from TRandom",20,0,20);
285 
286  double mu = 5;
287 
288  TF1 * f = new TF1("fps",poisson,1,0,1);
289  f->SetParameter(0,mu);
290 
292  TUnuran unr;
293 
294  // dari method (needs also the mode and pmf sum)
295  dist2.SetMode(int(mu) );
296  dist2.SetProbSum(1.0);
297  bool ret = unr.Init(dist2,"dari");
298  if (!ret) return;
299 
300  TStopwatch w;
301  w.Start();
302 
303  int n = NGEN;
304  for (int i = 0; i < n; ++i) {
305  int k = unr.SampleDiscr();
306  h1->Fill( double(k) );
307  }
308 
309  w.Stop();
310  cout << "Time using Unuran method " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
311 
312  w.Start();
313  for (int i = 0; i < n; ++i) {
314  h2->Fill( gRandom->Poisson(mu) );
315  }
316  cout << "Time using TRandom::Poisson " << "\t=\t\t " << w.CpuTime() << endl;
317 
318  c1->cd(++izone);
319  h1->SetMarkerStyle(20);
320  h1->Draw("E");
321  h2->Draw("same");
322 
323  std::cout << " chi2 test of UNURAN vs TRandom generated histograms: " << std::endl;
324  h1->Chi2Test(h2,"UUP");
325 
326 }
327 
328 //_____________________________________________
329 //
330 // example of empirical distributions
331 
332 void testEmpDistr() {
333 
334 
335  cout << "\nTest Empirical distributions using smoothing\n\n";
336 
337  // start with a set of data
338  // for example 1000 two-gaussian data
339  const int Ndata = 1000;
340  double x[Ndata];
341  for (int i = 0; i < Ndata; ++i) {
342  if (i < 0.5*Ndata )
343  x[i] = gRandom->Gaus(-1.,1.);
344  else
345  x[i] = gRandom->Gaus(1.,3.);
346  }
347 
348  TH1D * h0 = new TH1D("h0Ref","Starting data",100,-10,10);
349  TH1D * h1 = new TH1D("h1Unr","Unuran unbin Generated data",100,-10,10);
350  TH1D * h1b = new TH1D("h1bUnr","Unuran bin Generated data",100,-10,10);
351  TH1D * h2 = new TH1D("h2GR","Data from TH1::GetRandom",100,-10,10);
352 
353  h0->FillN(Ndata,x,0,1); // fill histogram with starting data
354 
355  TUnuran unr;
356  TUnuranEmpDist dist(x,x+Ndata,1);
357 
358 
359  TStopwatch w;
360  int n = NGEN;
361 
362  w.Start();
363  if (!unr.Init(dist)) return;
364  for (int i = 0; i < n; ++i) {
365  h1->Fill( unr.Sample() );
366  }
367 
368  w.Stop();
369  cout << "Time using Unuran unbin " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
370 
371  TUnuranEmpDist binDist(h0);
372 
373  w.Start();
374  if (!unr.Init(binDist)) return;
375  for (int i = 0; i < n; ++i) {
376  h1b->Fill( unr.Sample() );
377  }
378  w.Stop();
379  cout << "Time using Unuran bin " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
380 
381  w.Start();
382  for (int i = 0; i < n; ++i) {
383  h2->Fill( h0->GetRandom() );
384  }
385  cout << "Time using TH1::GetRandom " << "\t=\t\t " << w.CpuTime() << endl;
386 
387  c1->cd(++izone);
388 
389  h2->Draw();
390  h1->SetLineColor(kRed);
391  h1->Draw("same");
392  h1b->SetLineColor(kBlue);
393  h1b->Draw("same");
394 
395 
396 }
397 
398 
399 
400 void unuranDemo() {
401 
402  //gRandom->SetSeed(0);
403 
404  // load libraries
405  gSystem->Load("libMathCore");
406  gSystem->Load("libUnuran");
407 
408  // create canvas
409 
410  c1 = new TCanvas("c1_unuranMulti","Multidimensional distribution",10,10,1000,1000);
411  c1->Divide(2,4);
412  gStyle->SetOptFit();
413 
414  testStringAPI();
415  c1->Update();
416  testDistr1D();
417  c1->Update();
418  testDistrMultiDim();
419  c1->Update();
420  testDiscDistr();
421  c1->Update();
422  testEmpDistr();
423  c1->Update();
424 
425 
426 }
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
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
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:389
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:263
return c1
Definition: legend1.C:41
Definition: Rtypes.h:64
Random number generator class based on the maximally quidistributed combined Tausworthe generator b...
Definition: TRandom2.h:27
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
void SetMode(int mode)
set the mode of the distribution (location of maximum probability)
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1849
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
#define f(i)
Definition: RSha256.hxx:104
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition: TUnuran.cxx:382
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition: TUnuran.cxx:396
double sqrt(double)
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
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:1950
TUnuranDiscrDist class for one dimensional discrete distribution.
double breitwigner_pdf(double x, double gamma, double x0=0)
Probability density function of Breit-Wigner distribution, which is similar, just a different definit...
TH1F * h1
Definition: legend1.C:5
constexpr Double_t Pi()
Definition: TMath.h:38
TUnuranEmpDist class for describing empiral distributions.
double breitwigner_cdf(double x, double gamma, double x0=0)
Cumulative distribution function (lower tail) of the Breit_Wigner distribution and it is similar (jus...
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:304
A 3-Dim function with parameters.
Definition: TF3.h:28
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:1542
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH3.cxx:287
TUnuranMultiContDist class describing multi dimensional continuous distributions. ...
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
The Canvas class.
Definition: TCanvas.h:27
TUnuran class.
Definition: TUnuran.h:79
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
TUnuranContDist class describing one dimensional continuous distribution.
virtual void FillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride=1)
Fill this histogram with an array x and weights w.
Definition: TH1.cxx:3378
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:210
#define c(i)
Definition: RSha256.hxx:101
Definition: Rtypes.h:64
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:391
double exp(double)
virtual Double_t GetRandom() const
Return a random number distributed according the histogram bin contents.
Definition: TH1.cxx:4856
const Int_t n
Definition: legend1.C:16
Stopwatch class.
Definition: TStopwatch.h:28