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:559
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1855
#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:66
@ kBlue
Definition: Rtypes.h:66
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
R__EXTERN TStyle * gStyle
Definition: TStyle.h:414
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
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:649
A 3-Dim function with parameters.
Definition: TF3.h:28
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:617
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:2002
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3348
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3070
virtual Double_t GetRandom(TRandom *rng=nullptr) const
Return a random number distributed according the histogram bin contents.
Definition: TH1.cxx:4984
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:3451
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:312
Int_t Fill(Double_t) override
Invalid Fill method.
Definition: TH3.cxx:337
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 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 empirical distributions.
TUnuranMultiContDist class describing multi dimensional continuous distributions.
TUnuran class.
Definition: TUnuran.h:79
int SampleDiscr()
Sample discrete distributions.
Definition: TUnuran.cxx:422
const std::string & MethodName() const
used Unuran method
Definition: TUnuran.h:283
bool SampleMulti(double *x)
Sample multidimensional distributions.
Definition: TUnuran.cxx:436
bool Init(const std::string &distr, const std::string &method)
Initialize with Unuran string API interface.
Definition: TUnuran.cxx:75
double Sample()
Sample 1D distribution.
Definition: TUnuran.cxx:429
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...
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
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
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
constexpr Double_t Pi()
Definition: TMath.h:37
Author
Lorenzo Moneta

Definition in file unuranDemo.C.