Logo ROOT  
Reference Guide
unuranDemo.C File Reference

Detailed Description

Example macro to show unuran capabilities The results are compared with what is obtained using TRandom or TF1::GetRandom The macro is divided in 3 parts:

  • testStringAPI : show how to use string API of UNURAN to generate Gaussian random numbers
  • testDistr1D : show how to pass a 1D distribution object to UNURAN to generate numbers according to the given distribution object
  • testDistrMultiDIm : show how to pass a multidimensional distribution object to UNURAN

To execute the macro type in:

root[0]: gSystem->Load("libMathCore");
root[0]: gSystem->Load("libUnuran");
root[0]: .x unuranDemo.C+
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1850
#include "TStopwatch.h"
#include "TUnuran.h"
#include "TUnuranEmpDist.h"
#include "TH1.h"
#include "TH3.h"
#include "TF3.h"
#include "TMath.h"
#include "TRandom2.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "Math/ProbFunc.h"
#include "Math/DistFunc.h"
#include <iostream>
#include <cassert>
using std::cout;
using std::endl;
// number of distribution generated points
#define NGEN 1000000
int izone = 0;
TCanvas * c1 = 0;
// test using UNURAN string interface
void testStringAPI() {
TH1D * h1 = new TH1D("h1G","gaussian distribution from Unuran",100,-10,10);
TH1D * h2 = new TH1D("h2G","gaussian distribution from TRandom",100,-10,10);
cout << "\nTest using UNURAN string API \n\n";
TUnuran unr;
if (! unr.Init( "normal()", "method=arou") ) {
cout << "Error initializing unuran" << endl;
return;
}
int n = NGEN;
w.Start();
for (int i = 0; i < n; ++i) {
double x = unr.Sample();
h1->Fill( x );
}
w.Stop();
cout << "Time using Unuran method " << unr.MethodName() << "\t=\t " << w.CpuTime() << endl;
// use TRandom::Gaus
w.Start();
for (int i = 0; i < n; ++i) {
double x = gRandom->Gaus(0,1);
h2->Fill( x );
}
w.Stop();
cout << "Time using TRandom::Gaus \t=\t " << w.CpuTime() << endl;
assert(c1 != 0);
c1->cd(++izone);
h1->Draw();
c1->cd(++izone);
h2->Draw();
}
double distr(double *x, double *p) {
return ROOT::Math::breitwigner_pdf(x[0],p[0],p[1]);
}
double cdf(double *x, double *p) {
return ROOT::Math::breitwigner_cdf(x[0],p[0],p[1]);
}
// test of unuran passing as input a distribution object( a BreitWigner) distribution
void testDistr1D() {
cout << "\nTest 1D Continous distributions\n\n";
TH1D * h1 = new TH1D("h1BW","Breit-Wigner distribution from Unuran",100,-10,10);
TH1D * h2 = new TH1D("h2BW","Breit-Wigner distribution from GetRandom",100,-10,10);
TF1 * f = new TF1("distrFunc",distr,-10,10,2);
double par[2] = {1,0}; // values are gamma and mean
f->SetParameters(par);
TF1 * fc = new TF1("cdfFunc",cdf,-10,10,2);
fc->SetParameters(par);
// create Unuran 1D distribution object
// to use a different random number engine do:
TRandom2 * random = new TRandom2();
int logLevel = 2;
TUnuran unr(random,logLevel);
// select unuran method for generating the random numbers
std::string method = "tdr";
//std::string method = "method=auto";
// "method=hinv"
// set the cdf for some methods like hinv that requires it
// dist.SetCdf(fc);
//cout << "unuran method is " << method << endl;
if (!unr.Init(dist,method) ) {
cout << "Error initializing unuran" << endl;
return;
}
w.Start();
int n = NGEN;
for (int i = 0; i < n; ++i) {
double x = unr.Sample();
h1->Fill( x );
}
w.Stop();
cout << "Time using Unuran method " << unr.MethodName() << "\t=\t " << w.CpuTime() << endl;
w.Start();
for (int i = 0; i < n; ++i) {
double x = f->GetRandom();
h2->Fill( x );
}
w.Stop();
cout << "Time using TF1::GetRandom() \t=\t " << w.CpuTime() << endl;
c1->cd(++izone);
h1->Draw();
c1->cd(++izone);
h2->Draw();
std::cout << " chi2 test of UNURAN vs GetRandom generated histograms: " << std::endl;
h1->Chi2Test(h2,"UUP");
}
// 3D gaussian distribution
double gaus3d(double *x, double *p) {
double sigma_x = p[0];
double sigma_y = p[1];
double sigma_z = p[2];
double rho = p[2];
double u = x[0] / sigma_x ;
double v = x[1] / sigma_y ;
double w = x[2] / sigma_z ;
double c = 1 - rho*rho ;
double result = (1 / (2 * TMath::Pi() * sigma_x * sigma_y * sigma_z * sqrt(c)))
* exp (-(u * u - 2 * rho * u * v + v * v + w*w) / (2 * c));
return result;
}
// test of unuran passing as input a multi-dimension distribution object
void testDistrMultiDim() {
cout << "\nTest Multidimensional distributions\n\n";
TH3D * h1 = new TH3D("h13D","gaussian 3D distribution from Unuran",50,-10,10,50,-10,10,50,-10,10);
TH3D * h2 = new TH3D("h23D","gaussian 3D distribution from GetRandom",50,-10,10,50,-10,10,50,-10,10);
TF3 * f = new TF3("g3d",gaus3d,-10,10,-10,10,-10,10,3);
double par[3] = {2,2,0.5};
f->SetParameters(par);
//std::string method = "method=vnrou";
//std::string method = "method=hitro;use_boundingrectangle=false ";
std::string method = "hitro";
if ( ! unr.Init(dist,method) ) {
cout << "Error initializing unuran" << endl;
return;
}
w.Start();
double x[3];
for (int i = 0; i < NGEN; ++i) {
unr.SampleMulti(x);
h1->Fill(x[0],x[1],x[2]);
}
w.Stop();
cout << "Time using Unuran method " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
assert(c1 != 0);
c1->cd(++izone);
h1->Draw();
// need to set a reasonable number of points in TF1 to get acceptable quality from GetRandom to
int np = 200;
f->SetNpx(np);
f->SetNpy(np);
f->SetNpz(np);
w.Start();
for (int i = 0; i < NGEN; ++i) {
f->GetRandom3(x[0],x[1],x[2]);
h2->Fill(x[0],x[1],x[2]);
}
w.Stop();
cout << "Time using TF1::GetRandom \t\t=\t " << w.CpuTime() << endl;
c1->cd(++izone);
h2->Draw();
std::cout << " chi2 test of UNURAN vs GetRandom generated histograms: " << std::endl;
h1->Chi2Test(h2,"UUP");
}
//_____________________________________________
//
// example of discrete distributions
double poisson(double * x, double * p) {
return ROOT::Math::poisson_pdf(int(x[0]),p[0]);
}
void testDiscDistr() {
cout << "\nTest Discrete distributions\n\n";
TH1D * h1 = new TH1D("h1PS","Unuran Poisson prob",20,0,20);
TH1D * h2 = new TH1D("h2PS","Poisson dist from TRandom",20,0,20);
double mu = 5;
TF1 * f = new TF1("fps",poisson,1,0,1);
f->SetParameter(0,mu);
TUnuran unr;
// dari method (needs also the mode and pmf sum)
dist2.SetMode(int(mu) );
dist2.SetProbSum(1.0);
bool ret = unr.Init(dist2,"dari");
if (!ret) return;
w.Start();
int n = NGEN;
for (int i = 0; i < n; ++i) {
int k = unr.SampleDiscr();
h1->Fill( double(k) );
}
w.Stop();
cout << "Time using Unuran method " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
w.Start();
for (int i = 0; i < n; ++i) {
h2->Fill( gRandom->Poisson(mu) );
}
cout << "Time using TRandom::Poisson " << "\t=\t\t " << w.CpuTime() << endl;
c1->cd(++izone);
h1->Draw("E");
h2->Draw("same");
std::cout << " chi2 test of UNURAN vs TRandom generated histograms: " << std::endl;
h1->Chi2Test(h2,"UUP");
}
//_____________________________________________
//
// example of empirical distributions
void testEmpDistr() {
cout << "\nTest Empirical distributions using smoothing\n\n";
// start with a set of data
// for example 1000 two-gaussian data
const int Ndata = 1000;
double x[Ndata];
for (int i = 0; i < Ndata; ++i) {
if (i < 0.5*Ndata )
x[i] = gRandom->Gaus(-1.,1.);
else
x[i] = gRandom->Gaus(1.,3.);
}
TH1D * h0 = new TH1D("h0Ref","Starting data",100,-10,10);
TH1D * h1 = new TH1D("h1Unr","Unuran unbin Generated data",100,-10,10);
TH1D * h1b = new TH1D("h1bUnr","Unuran bin Generated data",100,-10,10);
TH1D * h2 = new TH1D("h2GR","Data from TH1::GetRandom",100,-10,10);
h0->FillN(Ndata,x,0,1); // fill histogram with starting data
TUnuran unr;
TUnuranEmpDist dist(x,x+Ndata,1);
int n = NGEN;
w.Start();
if (!unr.Init(dist)) return;
for (int i = 0; i < n; ++i) {
h1->Fill( unr.Sample() );
}
w.Stop();
cout << "Time using Unuran unbin " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
TUnuranEmpDist binDist(h0);
w.Start();
if (!unr.Init(binDist)) return;
for (int i = 0; i < n; ++i) {
h1b->Fill( unr.Sample() );
}
w.Stop();
cout << "Time using Unuran bin " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;
w.Start();
for (int i = 0; i < n; ++i) {
h2->Fill( h0->GetRandom() );
}
cout << "Time using TH1::GetRandom " << "\t=\t\t " << w.CpuTime() << endl;
c1->cd(++izone);
h2->Draw();
h1->Draw("same");
h1b->Draw("same");
}
void unuranDemo() {
//gRandom->SetSeed(0);
// load libraries
gSystem->Load("libMathCore");
gSystem->Load("libUnuran");
// create canvas
c1 = new TCanvas("c1_unuranMulti","Multidimensional distribution",10,10,1000,1000);
c1->Divide(2,4);
testStringAPI();
c1->Update();
testDistr1D();
c1->Update();
testDistrMultiDim();
c1->Update();
testDiscDistr();
c1->Update();
testEmpDistr();
c1->Update();
}
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
double sqrt(double)
double exp(double)
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
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:27
1-Dim function class
Definition: TF1.h:210
A 3-Dim function with parameters.
Definition: TF3.h:28
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual Double_t GetRandom() const
Return a random number distributed according the histogram bin contents.
Definition: TH1.cxx:4856
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
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
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
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:304
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:263
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:391
Stopwatch class.
Definition: TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
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 Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
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
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:382
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:396
bool Init(const std::string &distr, const std::string &method)
initialize with Unuran string interface
Definition: TUnuran.cxx:79
double Sample()
Sample 1D distribution User is responsible for having previously correctly initialized with TUnuran::...
Definition: TUnuran.cxx:389
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)
Definition: 3DDistances.cxx:48
constexpr Double_t Pi()
Definition: TMath.h:38
Author
Lorenzo Moneta

Definition in file unuranDemo.C.