Soluzioni: Lezione 4
Contents
Soluzioni: Lezione 4¶
4.1 il generatore lineare congruenziale¶
esercizio01.cpp
//c++ -o es1 esercizio01.cpp
#include <iostream>
static const int A = 214013;
static const int C = 2531011;
static const int M = 2147483647;
void rand(long int &num){
num = (A* num + C)%M;
return;
}
int main(){
long int seed_ = 0;
for(int i = 0; i < 10; ++i){
rand(seed_);
std::cout << i << ": " << seed_ << std::endl;
}
return 0;
}
4.2 il generatore lineare congruenziale implementato sotto forma di classe¶
esercizio02.cpp
//c++ -o es2 esercizio02.cpp
#include <iostream>
class LCG{
public:
LCG(): seed_(0) {};
LCG(int s): seed_(s) {};
~LCG() {};
void set_seed(long int s) { seed_ = s; };
long int* rnd();
private:
long int seed_; // unsigned long int if M > 2147483647
const int A = 214013;
const int C = 2531011;
const int M = 2147483647;
};
long int* LCG::rnd(){
seed_ = (A*seed_ + C)%M; //updating the seed x_{n+1} = (A*x_{n} + C)%M
return &seed_; //return reference to class member
}
int main(){
LCG rand; // initialize with seed = 0
for(int i = 0; i < 10; ++i){
std::cout << i << ": " << *rand.rnd() << std::endl;
}
return 0;
}
4.3 inizializzazione del seed di un generatore¶
esercizio04.cpp
//c++ -o es4 `root-config --glibs --cflags` esercizio04.cpp
#include <iostream>
#include "TH1F.h"
#include "TCanvas.h"
static const int A = 214013;
static const int C = 2531011;
static const int M = 2147483647;
long int seed_;
void rnd(){
/*
Updating the seed with LCG algorithm
*/
seed_ = (A* seed_ + C)%M;
return;
}
float rand_range(float min, float max){
/*
Updating the seed with LCG algorithm and scaling it in the required range [min, max]
*/
rnd();
return min + (max - min)*seed_ / static_cast<float> (M);
}
int main(){
seed_ = 0; //setting global seed of the sequence
float min = -3;
float max = 3;
TH1F histo ("histo", "prova", 5, min, max) ;
for(int i =0; i < 10; ++i){
float rand = rand_range(min, max);
std::cout << rand << std::endl;
histo.Fill (rand);
}
TCanvas c1 ;
histo.Draw () ;
c1.Print ("histo.png", "png") ;
return 0;
}
4.4 generazione di numeri pseudo-casuali uniformemente distribuiti¶
esercizio04.cpp
//c++ -o es4 `root-config --glibs --cflags` esercizio04.cpp
#include <iostream>
#include "TH1F.h"
#include "TCanvas.h"
static const int A = 214013;
static const int C = 2531011;
static const int M = 2147483647;
long int seed_;
void rnd(){
/*
Updating the seed with LCG algorithm
*/
seed_ = (A* seed_ + C)%M;
return;
}
float rand_range(float min, float max){
/*
Updating the seed with LCG algorithm and scaling it in the required range [min, max]
*/
rnd();
return min + (max - min)*seed_ / static_cast<float> (M);
}
int main(){
seed_ = 0; //setting global seed of the sequence
float min = -3;
float max = 3;
TH1F histo ("histo", "prova", 5, min, max) ;
for(int i =0; i < 10; ++i){
float rand = rand_range(min, max);
std::cout << rand << std::endl;
histo.Fill (rand);
}
TCanvas c1 ;
histo.Draw () ;
c1.Print ("histo.png", "png") ;
return 0;
}
4.5 generazione di numeri pseudo-casuali con il metodo try-and-catch¶
esercizio05.cpp
// c++ -o es5 `root-config --glibs --cflags` esercizio05.cpp
#include <iostream>
#include <cmath>
#include <iomanip>
#include "TH1F.h"
#include "TCanvas.h"
static const int A = 214013;
static const int C = 2531011;
static const int M = 2147483647;
long int seed_;
void rnd(){
/*
Updating the seed with LCG algorithm
*/
seed_ = (A* seed_ + C)%M;
return;
}
float rand_range(float min, float max){
/*
Updating the seed with LCG algorithm and scaling it in the required range [min, max]
*/
rnd();
return min + (max - min)*seed_ / static_cast<float> (M);
}
float fgaus (float x){
/*
Gaussian distribution sigma = 1
*/
return exp (-0.5 * x * x) ;
}
float fpois (float x){
/*
Possion distribution mu=2.
The definition of a factorial cannot be defined for non-integers,
but a common generalization is the gamma function, defined for all
positive reals and all negative non-integers -> true gamma in stl (std::tgamma).
*/
float fact = std::tgamma(x+1);
return exp(-3)*pow(3, x)/fact;
}
float funif (float x){
/*
Uniform distribution in -2, 2 -> constant at 1/4
*/
return float(1)/4;
}
float rand_TAC (float f (float), float xMin, float xMax, float yMax){
/*
Try and catch. input requires a function (f) taking as
inputs float values and output float values. The algorithm also
requires x_min and x_max to extract a uniform number in this range.
It then requires the y_max of the pdf (as pdf are non-negative y_min=0).
*/
float x = 0. ;
float y = 0. ;
do {
x = rand_range (xMin, xMax) ;
y = rand_range (0, yMax) ;
} while (y > f (x)) ;
return x ;
}
void visualize(float* vals, int dim, int bins, float x_min, float x_max){
/*
Simple visualization of the distribution in the shell output.
This function requires a pointer to the first element of an array of floats, its
dimension, the bin you want to represent and the edges x_min and x_max.
*/
float step = (x_max - x_min)/bins;
std::cout << std::left << std::setw(4) << "Bin Range" << std::right << std::setw(10) << " " << "Entries" << std::endl ;
for(int j = 0; j < bins; ++j){
float lower_edge = x_min + j*step;
float upper_edge = x_min + (j+1)*step;
std::cout << std::left << std::setw(4) << lower_edge << "," << std::right << std::setw(4) << upper_edge << ": " ;
for(int i = 0; i < dim; ++i){
if (lower_edge <= vals[i] && vals[i] < upper_edge) std::cout << "=";
}
std::cout << "\n";
}
return;
}
int main(){
seed_ = 0; //setting global seed of the sequence
int dim = 200; //extracting 200 values
float vals[dim]; //defining the array to same this pseudo-random values
float x_min = -3;
float x_max = 3;
TH1F hist_gauss ("hist_gauss", "GAUSS", 5, x_min, x_max) ;
for(int i=0; i < dim; ++i){
vals[i] = rand_TAC(fgaus, x_min, x_max, 10);
hist_gauss.Fill (vals[i]);
}
TCanvas c1 ;
hist_gauss.Draw () ;
c1.Print ("hist_gauss.png", "png") ;
//Rough visualization of gaussian
std::cout << "\n\nGAUSSIAN\n\n" << std::endl;
visualize(vals, dim, 10, x_min, x_max);
x_min = 0;
x_max = 10;
TH1F hist_poisson ("hist_poisson", "POISSON", 5, x_min, x_max) ;
for(int i=0; i < dim; ++i){
vals[i] = rand_TAC(fpois, x_min, x_max, 10);
hist_poisson.Fill (vals[i]);
}
TCanvas c2 ;
hist_poisson.Draw () ;
c2.Print ("hist_poisson.png", "png") ;
//Rough visualization of poisson
std::cout << "\n\nPOISSONIAN\n\n" << std::endl;
visualize(vals, dim, 10, x_min, x_max);
x_min = -2;
x_max = 2;
TH1F hist_uniform ("hist_uniform", "UNIFORM", 5, x_min, x_max) ;
for(int i=0; i < dim; ++i){
vals[i] = rand_TAC(funif, x_min, x_max, float(1)/4);
hist_uniform.Fill (vals[i]);
}
TCanvas c3 ;
hist_uniform.Draw () ;
c3.Print ("hist_uniform.png", "png") ;
//Rough visualization of uniform
std::cout << "\n\nUNIFORM\n\n" << std::endl;
visualize(vals, dim, 10, x_min, x_max);
return 0;
}
4.6 generazione di numeri pseudo-casuali con il metodo della funzione inversa¶
esercizio06.cpp
//c++ -o es6 `root-config --glibs --cflags` esercizio06.cpp
#include <iostream>
#include <cmath>
#include <iomanip>
#include "TH1F.h"
#include "TCanvas.h"
static const int A = 214013;
static const int C = 2531011;
static const int M = 2147483647;
long int seed_;
void rnd(){
/*
Updating the seed with LCG algorithm
*/
seed_ = (A* seed_ + C)%M;
return;
}
float rand_range(float min, float max){
/*
Updating the seed with LCG algorithm and scaling it in the required range [min, max]
*/
rnd();
return min + (max - min)*seed_ / static_cast<float> (M);
}
float inv_exp(float y){
/*
For completeness we define the \lambda parameter of the
exponential distribution (actually a waste of memory but cool to play).
pdf(x) = \lambda * exp(-\lambda x) x >= 0, 0 otherwise.
F(x) = int_{0}^{x} pdf(x)dx = 1 - exp(-\lambda * x) for x >= 0, 0 otherwise.
F^{-1}(y) = - (ln(1-y)) / \lambda
*/
float lambda = 1;
return - (log(1-y)/lambda);
}
void visualize(float* vals, int dim, int bins, float x_min, float x_max){
/*
Simple visualization of the distribution in the shell output.
This function requires a pointer to the first element of an array of floats, its
dimension, the bin you want to represent and the edges x_min and x_max.
*/
float step = (x_max - x_min)/bins;
std::cout << std::left << std::setw(4) << "Bin Range" << std::right << std::setw(10) << " " << "Entries" << std::endl ;
for(int j = 0; j < bins; ++j){
float lower_edge = x_min + j*step;
float upper_edge = x_min + (j+1)*step;
std::cout << std::left << std::setw(4) << lower_edge << "," << std::right << std::setw(4) << upper_edge << ": " ;
for(int i = 0; i < dim; ++i){
if (lower_edge <= vals[i] && vals[i] < upper_edge) std::cout << "=";
}
std::cout << "\n";
}
return;
}
int main(){
seed_ = 0; //setting global seed of the sequence
int dim = 200; //extracting 200 values
float vals[dim]; //defining the array to same this pseudo-random values
float y_min = 0;
float y_max = 1;
TH1F histo ("histo", "EXPONENTIAL", 10, 0, 5) ;
for(int i=0; i < dim; ++i){
vals[i] = inv_exp(rand_range(y_min, y_max));
histo.Fill(vals[i]);
}
//Rough visualization of gaussian
std::cout << "\n\nEXPONENTIAL\n\n" << std::endl;
visualize(vals, dim, 10, 0, 5);
TCanvas c1 ;
histo.Draw () ;
c1.Print ("EXPONENTIAL.png", "png") ;
}
4.7 generazione di numeri pseudo-casuali con il metodo del teorema centrale del limite¶
esercizio07.cpp
//c++ -o es7 `root-config --glibs --cflags` esercizio07.cpp
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <iomanip>
#include "TString.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TApplication.h"
static const int seed = 1234567;
float rand_range (float min, float max){
return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
}
float normal(int N, float min, float max){
/*
Basic principle:
Uniform distribution in [a,b] has a mean of 0.5*(a+b) and a variance 1/12 * (b-a)^2.
The mean of N uniformly distributed numbers is a variable that approximates a gaussian
distribution (exactly gaussian when N -> infinity).
*/
float mean = 0;
for(int i = 0; i < N; ++i) mean += rand_range(min, max);
return mean/N;
}
void visualize(float* vals, int dim, int bins, float x_min, float x_max){
/*
Simple visualization of the distribution in the shell output.
This function requires a pointer to the first element of an array of floats, its
dimension, the bin you want to represent and the edges x_min and x_max.
*/
float step = (x_max - x_min)/bins;
std::cout << std::left << std::setw(4) << "Bin Range" << std::right << std::setw(10) << " " << "Entries" << std::endl ;
for(int j = 0; j < bins; ++j){
float lower_edge = x_min + j*step;
float upper_edge = x_min + (j+1)*step;
std::cout << std::left << std::setw(4) << lower_edge << "," << std::right << std::setw(4) << upper_edge << ": " ;
for(int i = 0; i < dim; ++i){
if (lower_edge <= vals[i] && vals[i] < upper_edge) std::cout << "=";
}
std::cout << "\n";
}
return;
}
int main(){
srand(seed); //initializing seed for reproducibility of results
int Ns[4] = {1, 2, 10, 50};
int dim = 10000;
int nBins = 50;
float vals[dim];
// Computing the distributions
for(auto N : Ns){
float sum = 0;
float sum_q = 0;
for(int i = 0; i < dim; ++i){
vals[i] = normal(N, -3, 3);
sum += vals[i];
sum_q += pow(vals[i],2);
}
float mean = sum/dim;
float var = (sum_q + dim*pow(mean, 2) - 2*mean*sum)/(dim-1);
std::cout << "\n\nGAUSSIAN N,mean,var: " << N << " " << mean << " " << var << "\n\n" << std::endl;
//visualize(vals, dim, 10, -2, 2);
}
// Same thing but with standardization of the variables to obtain
// values distributed normally ~ N(0,1)
TApplication myApp("myApp", NULL, NULL);
TCanvas c0("c0","c0", 10, 10, 700, 500);
c0.Divide(2,2);
TH1F arrHisto[4];
for (unsigned int i = 0; i < 4; i++)
{
TString myStr = "GAUSSIAN STANDARD ";
myStr += Ns[i];
//std::cout << "Histo name : " << myStr << std::endl;
arrHisto[i] = TH1F(myStr, myStr, nBins, -3, 3);
}
for (unsigned int k = 0; k < 4; k++){
float sum = 0;
float sum_q = 0;
int N = Ns[k];
for(int i = 0; i < dim; ++i){
vals[i] = normal(N, -3, 3);
sum += vals[i];
sum_q += pow(vals[i],2);
}
float mean = sum/dim;
float var = (sum_q + dim*pow(mean, 2) - 2*mean*sum)/(dim-1);
//standardizing
// Z = (\hat{X} - \mu)/(\sigma) -> distributed normally ~ N(0,1)
sum = 0;
sum_q = 0;
for(int i = 0; i < dim; ++i){
vals[i] = (vals[i] - mean)/sqrt(var);
arrHisto[k].Fill(vals[i]);
sum += vals[i];
sum_q += pow(vals[i],2);
}
mean = sum/dim;
var = (sum_q + dim*pow(mean, 2) - 2*mean*sum)/(dim-1);
std::cout << "\n\nGAUSSIAN STANDARD N,mean,var: " << N << " " << mean << " " << var << "\n\n" << std::endl;
//visualize(vals, dim, 10, -3, 3);
}
for (unsigned int i = 0; i < 4; i++)
{
c0.cd(i+1);
arrHisto[i].Draw();
}
c0.Modified();
c0.Update();
myApp.Run();
return 0;
}
esercizio07_v2.cpp
/*
c++ -o es7_v2 esercizio07_v2.cpp `root-config --cflags --glibs`
*/
#include <iostream>
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TString.h"
#define NTOY 10000
#define MAXEXTR 20
#define NHISTO 20
double rand_range(const double min, const double max)
{
return min + (max - min) * (rand() / static_cast<double>(RAND_MAX));
}
int main()
{
srand(10);
TApplication myApp("myApp", NULL, NULL);
TCanvas c0("c0","c0", 10, 10, 700, 500);
c0.Divide(4,5);
// ############################
// # CONST INTERNAL VARIABLES #
// ############################
const double min = -10;
const double max = 10;
const int nBins = 100;
TH1D arrHisto[NHISTO];
for (unsigned int i = 0; i < NHISTO; i++)
{
TString myStr = "h";
myStr += i;
std::cout << "Histo name : " << myStr << std::endl;
arrHisto[i] = TH1D(myStr, myStr, nBins, -5, 5);
}
// #############
// # MEGA LOOP #
// #############
for (unsigned int k = 0; k < NHISTO; k++)
for (unsigned int i = 0; i < NTOY; i++)
{
double avg = 0;
double sum2 = 0;
unsigned int nExtr = (k+1) * MAXEXTR / NHISTO;
for (unsigned int j = 0; j < nExtr; j++)
{
double val = rand_range(min, max);
avg += val;
sum2 += val*val;
}
avg /= nExtr;
sum2 /= (nExtr * nExtr);
// Var[x] = E[x^2] - E^2[x]
// Var[x] = E[x^2]
// z = (avg - MEAN) / SIGMA
arrHisto[k].Fill(avg / sqrt(sum2));
}
// #################
// # VISUALIZATION #
// #################
for (unsigned int i = 0; i < NHISTO; i++)
{
c0.cd(i+1);
arrHisto[i].Draw();
}
c0.Modified();
c0.Update();
myApp.Run();
return 0;
}
4.8 implementazione e test di una classe `statistiche`¶
esercizio08_09.cpp
//c++ -o es_8_9 statistiche.cc esercizio08_09.cpp
#include <iostream>
#include <cmath>
#include "statistiche.h"
static const int A = 214013;
static const int C = 2531011;
static const int M = 2147483647;
long int seed_;
//----------------- Redefining random number generator (LCG) -----------------
void rnd(){
/*
Updating the seed with LCG algorithm
*/
seed_ = (A* seed_ + C)%M;
return;
}
float rand_range(float min, float max){
/*
Updating the seed with LCG algorithm and scaling it in the required range [min, max]
*/
rnd();
return min + (max - min)*seed_ / static_cast<float> (M);
}
//----------------- MAIN -----------------
int main(){
seed_ = 0; //setting global seed of the sequence
float x_min = -5;
float x_max = +5;
int dim = 100000;
statistiche* stats = new statistiche(dim/4);
for(int i = 0; i < dim; ++i) stats->aggiungiNumero(rand_range(x_min, x_max));
std::cout << "Uniform distribution in [" << x_min << "," << x_max << "] : " << std::endl;
std::cout << "Theoretical mean (0.5(a+b)): " << 0.5*(x_min + x_max)
<< " Computed mean: " << stats->media()
<< "+-" << stats->dev_standard_media(true)
<< std::endl;
std::cout << "Theoretical variance (1/12 (b-a)^2): " << (float(1)/12)*pow(x_max - x_min, 2)
<< " Computed variance: " << stats->varianza_1(true)
<< "+-" << stats->dev_standard_media_SQ(true)
<< std::endl;
//Checking for CLT
stats->reset();
int d = 20; //number of extractions before computing mean
for(int i = 0; i < dim; ++i){
double sum = 0;
for(int j = 0; j < d; ++j){
sum += rand_range(x_min, x_max);
}
stats->aggiungiNumero(sum/d);
}
std::cout << "\n\nCheck for CLT with N: " << d << std::endl;
std::cout << "Theoretical mean (0.5(a+b)): " << 0.5*(x_min + x_max)
<< " Computed mean: " << stats->media()
<< "+-" << stats->dev_standard_media(true)
<< std::endl;
std::cout << "Theoretical variance (1/12 (b-a)^2)/N: " << (float(1)/12)*pow(x_max - x_min, 2)/d
<< " Computed variance: " << stats->varianza_1(true)
<< "+-" << stats->dev_standard_media_SQ(true)
<< std::endl;
return 0;
}
statistiche.h
#ifndef statistiche_h
#define statistiche_h
#include <iostream>
#include <cmath>
//----------------- Statistic class -----------------
class statistiche
{
public:
statistiche (): dim(10), last_idx(0), arr (new float[dim]), sum (0.), sum_q (0.) {};
statistiche (int d): dim(d), last_idx(0), arr(dim ? new float[dim] : nullptr), sum (0.), sum_q (0.) {} ;
~statistiche () { delete[] arr; } ;
// aggiunge un numero alle informazioni salvate
void aggiungiNumero (float val);
// resituisce la media dei numeri aggiunti
double media () const ;
// restituisce la varianza dei numeri aggiunti
// in caso di opzione true, applica la correzione di Bessel
double varianza_1 (bool corretta = false) const ;
// restituisce la varianza dei numeri aggiunti
// in caso di opzione true, applica la correzione di Bessel
double varianza_2 (bool corretta = false) const ;
// resituisce la deviazione standard
// in caso di opzione true, applica la correzione di Bessel
double dev_standard (bool corretta = false) const ;
// resituisce la deviazione standard dalla media
double dev_standard_media (bool corretta = false) const ;
//restituisce varianza degli scarti quadratici
double varianza_SQ(bool corretta = false) const;
//restituisce std dev degli scarti quadratici
double dev_standard_SQ(bool corretta = false) const;
//restituisce std dev della media degli scarti quadratici
double dev_standard_media_SQ(bool corretta = false) const;
// resituisce il numero dei numeri aggiunti
double N () const { return last_idx; } ;
//reset
void reset() { for(int i = 0; i < last_idx; ++i) arr[i] = 0; last_idx = 0; sum = 0. ; sum_q = 0. ;};
private:
int dim;
int last_idx;
float* arr;
double sum ;
double sum_q ;
} ;
#endif
statistiche.cc
#include "statistiche.h"
void statistiche::aggiungiNumero (float val){
if (last_idx >= dim)
{
float * tempo = new float [2*dim] ;
for (int i = 0 ; i < dim ; ++i) tempo[i] = arr[i] ;
dim = dim * 2 ;
delete [] arr ;
arr = tempo ;
}
arr[last_idx] = val;
sum += val ;
sum_q += val * val ;
++last_idx;
return;
}
double statistiche::media() const{
/*
\hat{\mu} = sum(x_i)/N
*/
return sum/last_idx;
}
double statistiche::varianza_1(bool corretta) const{
/*
First variant. Based on the following.
Take the numerator of the variance formula:
sum[(x_i - \mu)^2] = sum((x_i)^2 + (\mu)^2 - 2x_i\mu) = \sum(x_i^2) + N\mu^2 - 2\mu sum(x_i)
two sums -> x_i and x_i^2. The first one is also useful to compute the mean while the second
is only instrumental for the variance computation.
*/
if (last_idx < 2) return 0 ;
double mean = media () ;
if (corretta) return (sum_q + last_idx * mean * mean - 2 * mean * sum) / (last_idx - 1) ;
else return sum_q / last_idx - mean * mean ;
}
double statistiche::varianza_2(bool corretta) const{
/*
Second variant. Performs two scans of the array -> worse performances.
First computes the mean then sums (x_i - \mu)^2. Lastly takes the mean
of the square deviations (with or without Bessel correction).
*/
if (last_idx < 2) return 0 ;
double mean = media();
double sum_q = 0;
for(int i = 0; i < last_idx; ++i){
sum_q += pow(arr[i] - mean, 2);
}
if (corretta) return sum_q/(last_idx-1);
else return sum_q/last_idx;
}
double statistiche::dev_standard (bool corretta) const{
/*
dev_std = sqrt(var)
*/
return sqrt(varianza_1(corretta));
}
double statistiche::dev_standard_media (bool corretta) const{
/*
dev_std_media = (dev_std)/sqrt(N)
*/
return sqrt(varianza_1(corretta)/last_idx);
}
double statistiche::varianza_SQ(bool corretta) const{
double var = varianza_1(); // mean of the squared deviations var = (x_i - \mu)^2 / N
double mean = media(); // \mu
if (corretta) return sum_q/(last_idx - 1);
else return sum_q/last_idx;
}
double statistiche::dev_standard_SQ (bool corretta) const{
/*
dev_std = sqrt(var)
*/
return sqrt(varianza_SQ(corretta));
}
double statistiche::dev_standard_media_SQ (bool corretta) const{
/*
dev_std_media = sqrt(var)
*/
return sqrt(varianza_SQ(corretta)/last_idx);
}