Soluzioni: Lezione 10
Contents
Soluzioni: Lezione 10¶
10.1 - 10.3 determinazione del massimo della verosimiglianza¶
esercizio_01.cpp
/*
c++ -o esercizio_01 funzioni.cc casual.cc esercizio_01.cpp
Testo:
Si scriva una libreria di funzioni per determinare il parametro τ della distribuzione
esponenziale utilizzata per generare gli eventi presenti nel file dati_esponenziali.txt,
a partire dal vector di numeri creato negli esercizi precedenti.
Si implementino i prototipi delle funzioni presentate a lezione.
Si confronti il risultato ottenuto con la media dei numeri salvati nel vector.
Come dipende il risultato dall'intervallo inziale passato alla funzione sezione_aurea_max?
*/
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include "statistiche_vector.h"
#include "funzioni.h"
#include "casual.h"
using namespace std ;
int main (int argc, char ** argv)
{
if (argc < 2)
{
cerr << "uso: " << argv[0] << " numero_di_eventi [t_zero]" << endl ;
exit (1) ;
}
double t_zero = 5. ;
if (argc == 3)
{
t_zero = atof (argv[2]) ;
}
vector<double> data ;
for (int i = 0 ; i < atoi (argv[1]) ; ++i)
{
data.push_back (rand_exp (t_zero)) ;
}
cout << "generati " << data.size () << " eventi" << endl ;
double media_v = media (data) ;
cout << "media = " << media_v << endl ;
double massimo = sezione_aurea_max (loglikelihood, 0.5 * media_v, 1.5 * media_v, data) ;
cout << "massimo della likelihood = " << massimo << endl ;
return 0 ;
}
casual.h
#ifndef casual_h
#define casual_h
double rand_range (double min, double max) ;
float rand_exp (double t_zero) ;
float rand_pois (double media) ;
double rand_IFC_Exp (double mu) ;
#endif
casual.cc
#include "casual.h"
#include <cstdlib>
#include <cmath>
double rand_range (double min, double max)
{
return min + (max - min) * rand () / static_cast<double> (RAND_MAX);
}
float rand_exp (double t_zero)
{
return -1. * log (1 - rand_range (0., 1.)) * t_zero ;
}
float rand_pois (double media)
{
double t_tot = rand_exp (1.) ;
int N_evt = 0 ;
while (t_tot < media)
{
++N_evt ;
t_tot += rand_exp (1.) ;
}
return N_evt ;
}
// generazione numeri casuali con il metodo dell'inversa della funzione cumulativa
double rand_IFC_Exp (double mu)
{
double y = rand_range (0., 1.) ;
double x = -mu * log (1 - y) ;
return x ;
}
funzioni.h
#ifndef funzioni_h
#define funzioni_h
#include <vector>
double esponenziale (double x, double tau) ;
double loglikelihood (
const std::vector<double> & data,
double param
) ;
double loglikelihoodprod (
const std::vector<double> & data,
double param
) ;
double sezione_aurea_max (
double logl (const std::vector<double> & , double),
double xMin,
double xMax,
const std::vector<double> & data,
double precision = 0.0001
) ;
double h (
const std::vector<double> & data,
double param,
double max
) ;
double bisezione (
double h (const std::vector<double> & , double, double),
double xMin,
double xMax,
const std::vector<double> & data,
double massimo,
double precision = 0.0001
) ;
#endif
funzioni.cc
#include "funzioni.h"
#include <cmath>
#include <iostream>
using namespace std ;
double esponenziale (double x, double tau)
{
if (tau == 0.) return 1. ;
return exp (-1. * x / tau) / tau ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double loglikelihood (
const vector<double> & data,
double param
)
{
double result = 0. ;
for (int i = 0 ; i < data.size () ; ++i) result += log (esponenziale (data.at (i), param)) ;
return result ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double loglikelihoodprod (
const vector<double> & data,
double param
)
{
double result = 1. ;
for (int i = 0 ; i < data.size () ; ++i) result *= esponenziale (data.at (i), param) ;
return log (result) ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double sezione_aurea_max (
double logl (const vector<double> & , double),
double xMin,
double xMax,
const vector<double> & data,
double precision
)
{
double x_2 = xMin ;
double x_3 = xMax ;
while ((xMax - xMin) > precision)
{
x_2 = xMin + (xMax - xMin) * 0.618 ;
x_3 = xMin + (xMax - xMin) * 0.382 ;
if (logl (data, x_3) < logl (data, x_2)) xMin = x_3 ;
else xMax = x_2 ;
// cout << "estremi " << xMin << "\t" << xMax << "\t" << xMax - xMin << "\n" ;
}
return 0.5 * (xMax + xMin) ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double h (
const vector<double> & data,
double param,
double max
)
{
return loglikelihood (data, param) + 0.5 - loglikelihood (data, max) ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double bisezione (
double h (const vector<double> & , double, double),
double xMin,
double xMax,
const vector<double> & data,
double massimo,
double precision
)
{
double xAve = xMin ;
while ((xMax - xMin) > precision)
{
xAve = 0.5 * (xMax + xMin) ;
if (h (data, xAve, massimo) * h (data, xMin, massimo) > 0.) xMin = xAve ;
else xMax = xAve ;
}
return xAve ;
}
10.4 varianza dello stimatore ottenuto¶
esercizio_04.cpp
/*
c++ -o esercizio_04 funzioni.cc casual.cc esercizio_04.cpp
Testo:
Si utilizzi il metodo della bisezione per trovare i due punti τ - στ e τ + στ nel caso
dell'esercizio precedente.
Si confronti il valore di στ ottenuto in questo modo con quello calcolato a partire
dalla media aritmetica dei numeri salvati nel vector
*/
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include "statistiche_vector.h"
#include "funzioni.h"
#include "casual.h"
using namespace std ;
int main (int argc, char ** argv)
{
if (argc < 2)
{
cerr << "uso: " << argv[0] << " numero_di_eventi [t_zero]" << endl ;
exit (1) ;
}
double t_zero = 5. ;
if (argc == 3)
{
t_zero = atof (argv[2]) ;
}
vector<double> data ;
for (int i = 0 ; i < atoi (argv[1]) ; ++i)
{
data.push_back (rand_exp (t_zero)) ;
}
cout << "generati " << data.size () << " eventi" << endl ;
double media_v = media (data) ;
cout << "media = " << media_v << endl ;
double massimo = sezione_aurea_max (loglikelihood, 0.5 * media_v, 1.5 * media_v, data) ;
cout << "massimo della likelihood = " << massimo << endl ;
double zero_sx = bisezione (h, 0.5 * media_v, massimo, data, massimo) ;
double zero_dx = bisezione (h, massimo, 1.5 * media_v, data, massimo) ;
cout << "zero_sx = " << zero_sx << endl ;
cout << "zero_dx = " << zero_dx << endl ;
cout << "sigma = " << 0.5 * (zero_dx - zero_sx) << endl ;
cout << "sigma stimata = " << media_v / sqrt (data.size ()) << endl ;
return 0 ;
}
statistiche_vector.h
#ifndef statistiche_vector_h
#define statistiche_vector_h
#include <vector>
template <class T>
T media (const std::vector<T> & input_v)
{
T somma = 0 ;
for (int i = 0 ; i < input_v.size () ; ++i) somma += input_v.at (i) ;
return somma / static_cast<float> (input_v.size ()) ;
}
template <class T>
T varianza (const std::vector<T> & input_v)
{
T somma = 0 ;
T sommaSq = 0 ;
for (int i = 0 ; i < input_v.size () ; ++i)
{
somma += input_v.at (i) ;
sommaSq += input_v.at (i) * input_v.at (i) ;
}
return sommaSq / static_cast<float> (input_v.size ()) -
(somma / static_cast<float> (input_v.size ()) * somma / static_cast<float> (input_v.size ())) ;
}
#endif
10.5 distribuzione di probabilità dello stimatore ottenuto¶
esercizio_05.cpp
/*
c++ -o esercizio_05 `root-config --glibs --cflags` funzioni.cc casual.cc esercizio_05.cpp
Testo:
Utilizzando il generatore di numeri pseudo-casuali secondo una pdf esponenziale sviluppato nella
Lezione 4, si disegni la distribuzione di probabilita' dello stimatore di τ in funzione del numero
di eventi a disposizione per la stima.
*/
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include "TH1F.h"
#include "TCanvas.h"
#include "TApplication.h"
#include "statistiche_vector.h"
#include "funzioni.h"
#include "casual.h"
using namespace std ;
int main (int argc, char ** argv)
{
if (argc < 2)
{
cerr << "uso: " << argv[0] << " mu_true numero_di_eventi [numero_di_toy = 10000]" << endl ;
exit (1) ;
}
double mu_true = atof (argv[1]) ;
int numero_eventi = atoi (argv[2]) ;
int N_toys = 10000 ;
if (argc > 3) N_toys = atoi (argv[3]) ;
TApplication theapp ("theapp", 0, 0) ;
TH1F h_max (
"h_max", "max likelihood distribution",
N_toys / 100,
mu_true - 3 * mu_true / sqrt (numero_eventi),
mu_true + 3 * mu_true / sqrt (numero_eventi)
) ;
for (int i_toy = 0 ; i_toy < N_toys ; ++i_toy)
{
if (i_toy % 1000 == 0) cout << "running toy " << i_toy << endl ;
vector<double> data_loc ;
for (int i_sample = 0 ; i_sample < numero_eventi ; ++i_sample)
{
data_loc.push_back (rand_IFC_Exp (mu_true)) ;
}
double media_v = media (data_loc) ;
double sigma_subsample = media_v / sqrt (data_loc.size ()) ;
double massimo = sezione_aurea_max (loglikelihood, 0.5 * media_v, 1.5 * media_v, data_loc) ;
h_max.Fill (massimo) ;
}
h_max.SetLineColor (kGray) ;
h_max.SetFillColor (kOrange + 1) ;
h_max.GetXaxis ()->SetTitle ("stimatore di #tau") ;
h_max.GetYaxis ()->SetTitle ("conteggi per bin") ;
TCanvas c1 ;
c1.SetLeftMargin (0.10) ;
h_max.Draw () ;
c1.Update();
theapp.Run () ;
return 0 ;
}
10.6 limite asintotico¶
esercizio_06.cpp
/*
c++ -o esercizio_06 `root-config --glibs --cflags` funzioni.cc casual.cc esercizio_06.cpp
Testo:
In regime asintotico, la distribuzione degli scarti (τ - τvero) / στ ha forma Normale.
Si utilizzi il metodo dei toy experiment per riempire l'istogramma degli scarti, dato un numero di eventi per toy experiment.
Si calcoli la media e la sigma della distribuzione degli scarti e se ne disegni il valore al variare del numero di eventi a
disposizione per la stima, riempiendo un TGraph con il numero di eventi a disposizione sull'asse orizziontale ed il valore
del parametro sull'asse verticale.
*/
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include "TH1F.h"
#include "TCanvas.h"
#include "TApplication.h"
#include "TGraph.h"
#include "statistiche_vector.h"
#include "funzioni.h"
#include "casual.h"
using namespace std ;
int main (int argc, char ** argv)
{
if (argc < 2)
{
cerr << "uso: " << argv[0] << " mu_true numero_di_eventi [numero_di_estrazioni = 10000]" << endl ;
exit (1) ;
}
double mu_true = atof (argv[1]) ;
int numero_eventi = 10000;
if (argc > 2) numero_eventi = atoi (argv[2]) ;
std::vector<int> N_toys_vec = {50, 100, 1000, 10000} ;
TApplication theapp ("theapp", 0, 0) ;
TGraph g;
g.SetMarkerStyle(8);
g.SetLineStyle(9);
g.SetTitle(";Ntoys;(#tau - #hat{#tau})/#hat{#sigma_{#tau}}");
std::vector<TH1F> histos;
for(int idx = 0; idx < N_toys_vec.size(); ++idx)
{
int N_toys = N_toys_vec.at(idx);
std::string h_name = "N_toy " + std::to_string(N_toys);
TH1F h(
h_name.c_str(), h_name.c_str(),
sqrt(N_toys),
-2,
2
) ;
h.SetLineColor(kBlue);
h.SetFillColor(0);
h.SetLineWidth(2);
h.GetXaxis()->SetTitle("(#hat{#tau} - #tau) / #sigma_{#tau}");
h.GetYaxis()->SetTitle("Events per bin");
std::vector<double> pulls;
for (int i_toy = 0 ; i_toy < N_toys ; ++i_toy)
{
if (i_toy % 1000 == 0) cout << "running toy " << i_toy << endl ;
vector<double> data_loc ;
for (int i_sample = 0 ; i_sample < numero_eventi ; ++i_sample)
{
data_loc.push_back (rand_IFC_Exp (mu_true)) ;
}
double media_v = media (data_loc) ;
double sigma_subsample = media_v / sqrt (data_loc.size ()) ;
double pull = (mu_true - media_v) / sigma_subsample;
pulls.push_back(pull);
h.Fill(pull);
}
g.SetPoint(idx, N_toys, media(pulls));
histos.push_back(h);
}
TCanvas c1 ;
c1.Divide(2,2 , 0.0001, 0.0001);
c1.SetLeftMargin (0.10) ;
for(int i = 0; i < N_toys_vec.size(); ++i){
c1.cd(i+1);
histos.at(i).Draw("hist");
};
c1.Update();
TCanvas c2 ;
c2.SetLogx();
c2.SetLeftMargin (0.10) ;
g.Draw("APL");
c2.Update();
theapp.Run () ;
return 0 ;
}