Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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"
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
51using std::cout;
52using std::endl;
53
54// number of distribution generated points
55#define NGEN 1000000
56
57int izone = 0;
58TCanvas * c1 = 0;
59
60// test using UNURAN string interface
61void 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
108double distr(double *x, double *p) {
109 return ROOT::Math::breitwigner_pdf(x[0],p[0],p[1]);
110}
111
112double 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
117void 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
189double 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
205void 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
275double poisson(double * x, double * p) {
276 return ROOT::Math::poisson_pdf(int(x[0]),p[0]);
277}
278
279void 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
332void 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();
391 h1->Draw("same");
392 h1b->SetLineColor(kBlue);
393 h1b->Draw("same");
394
395
396}
397
398
399
400void 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}
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
double sqrt(double)
double exp(double)
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
R__EXTERN TStyle * gStyle
Definition TStyle.h:412
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
static struct mg_connection * fc(struct mg_context *ctx)
Definition civetweb.c:3728
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
The Canvas class.
Definition TCanvas.h:23
1-Dim function class
Definition TF1.h:213
A 3-Dim function with parameters.
Definition TF3.h:28
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:618
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:2005
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
virtual Double_t GetRandom(TRandom *rng=nullptr) const
Return a random number distributed according the histogram bin contents.
Definition TH1.cxx:4942
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:3453
3-D histogram with a double per channel (see TH1 documentation)}
Definition TH3.h:305
Int_t Fill(Double_t)
Invalid Fill method.
Definition TH3.cxx:287
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition TRandom2.h:27
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:274
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:402
Stopwatch class.
Definition TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
void Stop()
Stop the stopwatch.
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:1541
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition TSystem.cxx:1853
TUnuranContDist class describing one dimensional continuous distribution.
TUnuranDiscrDist class for one dimensional discrete distribution.
void SetMode(int mode)
set the mode of the distribution (location of maximum probability)
void SetProbSum(double sum)
set the value of the sum of the probabilities in the given domain
TUnuranEmpDist class for describing empiral distributions.
TUnuranMultiContDist class describing multi dimensional continuous distributions.
TUnuran class.
Definition TUnuran.h:79
int SampleDiscr()
Sample discrete distributions User is responsible for having previously correctly initialized with TU...
Definition TUnuran.cxx:378
const std::string & MethodName() const
used Unuran method
Definition TUnuran.h:237
bool SampleMulti(double *x)
Sample multidimensional distributions User is responsible for having previously correctly initialized...
Definition TUnuran.cxx:392
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition TUnuran.cxx:75
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition TUnuran.cxx:385
double breitwigner_pdf(double x, double gamma, double x0=0)
Probability density function of Breit-Wigner distribution, which is similar, just a different definit...
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
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...
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5
double dist(Rotation3D const &r1, Rotation3D const &r2)
constexpr Double_t Pi()
Definition TMath.h:37