8.5. Solutions: Lecture 8#
8.5.1. Exercise 8.1#
es_8.1.py
#!/usr/bin/python
'''
Si scriva un programma che, fissato un numero *NMAX*,
generi *NToys* toy experiment,
ciascuno contenente un campione di *NMAX* eventi che seguano una distribuzione a scelta
e ne calcoli la media.
python3 es_8.1.py 100 1000
'''
import sys
from myrand import generate_uniform
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 3 :
print ('usage: ', sys.argv[0], 'N_toy N_evt')
exit ()
N_toy = int (argv[1])
N_evt = int (argv[2])
# loop over toys
for i in range (N_toy):
randlist = generate_uniform (N_evt)
toy_stats = stats (randlist)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.2. Exercise 8.2#
es_8.2.py
#!/usr/bin/python
'''
Si aggiunga al programma precedente un istogramma
che visualizzi la distribuzione delle medie al variare dei toy experiment.
python3 es_8.2.py 100 1000
'''
import sys
import numpy as np
from math import floor
import matplotlib.pyplot as plt
from myrand import generate_uniform
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 3 :
print ('usage: ', sys.argv[0], 'N_toy N_evt')
exit ()
N_toy = int (sys.argv[1])
N_evt = int (sys.argv[2])
means = []
# loop over toys
for i in range (N_toy):
randlist = generate_uniform (N_evt)
toy_stats = stats (randlist)
means.append (toy_stats.mean ())
means_stats = stats (means)
xMin = means_stats.mean () - 5. * means_stats.sigma ()
xMax = means_stats.mean () + 5. * means_stats.sigma ()
nBins = floor (len (means) / 10.) + 1 # number of bins of the histogram
bin_edges = np.linspace (xMin, xMax, nBins + 1) # edges o the histogram bins
fig, ax = plt.subplots ()
ax.set_title ('Histogram of the sigma of the mean over ' + str (N_toy) + ' toys', size=14)
ax.set_xlabel ('mean value')
ax.set_ylabel ('toys in bin')
ax.hist (means, # list of numbers
bins = bin_edges,
color = 'orange',
# normed = True,
)
plt.savefig ('es_8.2.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.3. Exercise 8.3#
es_8.3.cpp
#!/usr/bin/python
'''
Si utilizzi la classe ```stats``` sviluppata nelle lezioni precedenti
per confrontare la deviazione standard della media calcolata per ogni singolo toy
con la deviazione standard del campione delle medie.
python3 es_8.3.py 100 1000
'''
import sys
import numpy as np
from math import floor
import matplotlib.pyplot as plt
from myrand import generate_uniform
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 3 :
print ('usage: ', sys.argv[0], 'N_toy N_evt')
exit ()
N_toy = int (sys.argv[1])
N_evt = int (sys.argv[2])
means = []
sigma_means = []
# loop over toys
for i in range (N_toy):
randlist = generate_uniform (N_evt)
toy_stats = stats (randlist)
means.append (toy_stats.mean ())
sigma_means.append (toy_stats.sigma_mean ())
# compare the distribution of the sigma on the mean
# calculated for each toy to the sigma of the mean distribution
sigma_means_stats = stats (sigma_means)
means_stats = stats (means)
print ('sigma of the means disitribution: ', means_stats.sigma ())
print ('mean of the sigma of the means disitribution: ', sigma_means_stats.mean ())
# plot the distribution of the sigma on the mean
# calculated for each toy
xMin = sigma_means_stats.mean () - 5. * sigma_means_stats.sigma ()
xMax = sigma_means_stats.mean () + 5. * sigma_means_stats.sigma ()
nBins = floor (len (sigma_means) / 10.) + 1 # number of bins of the histogram
bin_edges = np.linspace (xMin, xMax, nBins + 1) # edges o the histogram bins
fig, ax = plt.subplots ()
ax.set_title ('Histogram of the sigma on the mean over ' + str (N_toy) + ' toys', size=14)
ax.set_xlabel ('mean value')
ax.set_ylabel ('toys in bin')
ax.hist (sigma_means, # list of numbers
bins = bin_edges,
color = 'orange',
# normed = True,
)
# get the sigma of the means distribution
ax.plot ([means_stats.sigma (), means_stats.sigma ()], plt.ylim ())
plt.savefig ('es_8.3.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.4. Exercise 8.4#
es_8.4.py
#!/usr/bin/python
'''
Si utilizzino due scatter plot per confrontare l'evoluzione
della deviazione standard della media calcolata per ogni singolo toy
con la deviazione standard del campione delle medie
al variare del numero di eventi generati in un singolo toy experiment.
python3 es_8.4.py 100 100 10000
'''
import sys
import numpy as np
from math import floor
import matplotlib.pyplot as plt
from myrand import generate_uniform
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 4 :
print ('usage: ', sys.argv[0], 'N_toy N_evt_min N_evt_max')
exit ()
N_toys = int (sys.argv[1])
N_evt_min = int (sys.argv[2])
N_evt_max = int (sys.argv[3])
# deviazione standard del campione delle medie
means_sigma = []
# single sigma of the mean
sigma_of_the_mean = []
sigmas_of_the_mean = []
N_events = []
X = []
while (N_evt_min < N_evt_max) :
means = []
# loop over toys
for i in range (N_toys):
randlist = generate_uniform (N_evt_min)
toy_stats = stats (randlist)
means.append (toy_stats.mean ())
if i == 0 : sigma_of_the_mean.append (toy_stats.sigma_mean ())
sigmas_of_the_mean.append (toy_stats.sigma_mean ())
X += [N_evt_min]
N_events.append (N_evt_min)
toy_stats = stats (means)
means_sigma.append (toy_stats.sigma ())
N_evt_min = N_evt_min * 2
fig, ax = plt.subplots ()
ax.set_title ('sigma of the means', size=14)
ax.set_xlabel ('number of events')
ax.set_ylabel ('sigma of the mean')
ax.scatter(X,sigmas_of_the_mean,label='stddev of the mean')
ax.plot (N_events, means_sigma, color = 'red', label = 'all toys')
ax.plot (N_events, sigma_of_the_mean, color = 'blue', label = 'single toy')
ax.set_xscale ('log')
ax.legend ()
plt.savefig ('es_8.4.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.5. Exercise 8.5#
es_8.5.py
#!/usr/bin/python
'''
Si implementi il metodo di integrazione hit-or-miss
con la funzione di esempio *f(x) = sin (x)*.
* Si scriva l'algoritmo che calcola l'integrale come una funzione esterna al programma ```main```,
facendo in modo che prenda come parametri di ingresso,
oltre agli estremi lungo l'asse *x* e l'asse *y*,
anche il numero di punti pseudo-casuali da generare.
* Si faccia in modo che l'algoritmo ritorni due elementi:
il primo elemento sia il valore dell'integrale,
il secondo sia la sua incertezza.
python3 es_8.5.py 0. 6.28 2. 1000
'''
import sys
import numpy as np
from integral import integral_HOM
def func (x) :
return 1. + np.sin (x) ;
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 5 :
print ('usage: ', sys.argv[0], 'xMin xMax yMin yMax N_evt')
exit ()
xMin = float (sys.argv[1])
xMax = float (sys.argv[2])
yMax = float (sys.argv[3])
N_evt = int (sys.argv[4])
integral, integral_unc = integral_HOM (func, xMin, xMax, yMax, N_evt)
print ('integral: ', integral)
print ('integral uncertainty: ', integral_unc)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
integral.py
#!/usr/bin/python
from myrand import generate_range, rand_range
from math import sqrt
def integral_HOM (func, xMin, xMax, yMax, N_evt) :
x_coord = generate_range (xMin, xMax, N_evt)
y_coord = generate_range (0., yMax, N_evt)
points_under = 0
for x, y in zip (x_coord, y_coord):
if (func (x) > y) : points_under = points_under + 1
A_rett = (xMax - xMin) * yMax
frac = float (points_under) / float (N_evt)
integral = A_rett * frac
integral_unc = A_rett**2 * frac * (1 - frac) / N_evt
return integral, integral_unc
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def integral_CrudeMC (g, xMin, xMax, N_rand) :
somma = 0.
sommaQ = 0.
for i in range (N_rand) :
x = rand_range (xMin, xMax)
somma += g(x)
sommaQ += g(x) * g(x)
media = somma / float (N_rand)
varianza = sommaQ /float (N_rand) - media * media
varianza = varianza * (N_rand - 1) / N_rand
lunghezza = (xMax - xMin)
return media * lunghezza, sqrt (varianza / float (N_rand)) * lunghezza
8.5.6. Exercise 8.6#
es_8.6.py
#!/usr/bin/python
'''
Si inserisca il calcolo dell'integrale dell'esercizio precedente in un ciclo che,
al variare del numero *N* di punti generati, mostri il valore dell'integrale
e della sua incertezza.
* Si utilizzi uno scatter plot per disegnare gli andamenti del valore dell'integrale
e della sua incertezza, al variare di *N* con ragione logaritmica.
python3 es_8.6.py 0. 6.28 2. 100 10000
'''
import sys
import numpy as np
import matplotlib.pyplot as plt
from integral import integral_HOM
def func (x) :
return 1. + np.sin (x) ;
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 6 :
print ('usage: ', sys.argv[0], 'xMin xMax yMin yMax N_evt_min N_evt_max')
exit ()
xMin = float (sys.argv[1])
xMax = float (sys.argv[2])
yMax = float (sys.argv[3])
N_evt_min = int (sys.argv[4])
N_evt_max = int (sys.argv[5])
N_events = []
integrals = []
integral_errors = []
while (N_evt_min < N_evt_max) :
integral, integral_unc = integral_HOM (func, xMin, xMax, yMax, N_evt_min)
integrals.append (integral)
integral_errors.append (integral_unc)
N_events.append (N_evt_min)
N_evt_min = N_evt_min * 2
fig, ax = plt.subplots ()
ax.set_title ('integral precision', size=14)
ax.set_xlabel ('number of events')
ax.set_ylabel ('integral and its error')
ax.errorbar (N_events, integrals, xerr = 0.0, yerr = integral_errors)
plt.savefig ('es_8.6.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.7. Exercise 8.7#
es_8.7.py
#!/usr/bin/python
'''
Si implementi il metodo di integrazione crude-MC
con la funzione di esempio *f(x) = sin (x)*.
* Si scriva l'algoritmo che calcola l'integrale come una funzione esterna al programma ```main```,
facendo in modo che prenda come parametri di ingresso,
oltre agli estremi lungo l'asse *x*,
anche il numero di punti pseudo-casuali da generare.
* Si faccia in modo che l'algoritmo ritorni un contenitore contenente due elementi:
il primo elemento sia il valore dell'integrale,
il secondo sia la sua incertezza.
python3 es_8.7.py 0. 6.28 2. 1000
'''
import sys
import numpy as np
from integral import integral_CrudeMC
def func (x) :
return 1. + np.sin (x) ;
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 5 :
print ('usage: ', sys.argv[0], 'xMin xMax yMin yMax N_evt')
exit ()
xMin = float (sys.argv[1])
xMax = float (sys.argv[2])
yMax = float (sys.argv[3])
N_evt = int (sys.argv[4])
integral, integral_unc = integral_CrudeMC (func, xMin, xMax, N_evt)
print ('integral: ', integral)
print ('integral uncertainty: ', integral_unc)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.8. Exercise 8.8#
es_8.8.py
#!/usr/bin/python
'''
Si inserisca il calcolo dell'integrale dell'esercizio precedente in un ciclo che,
al variare del numero *N* di punti generati, mostri il valore dell'integrale
e della sua incertezza.
* Si disegnino gli andamenti del valore dell'integrale
e della sua incertezza, al variare di *N* con ragione logaritmica.
* Si sovrapponga questo andamento a quello ottenuto dallo svolgimento dell'Esercizio 8.6.
python3 es_8.8.py 0. 6.28 2. 100 10000
'''
import sys
import numpy as np
import matplotlib.pyplot as plt
from integral import integral_HOM, integral_CrudeMC
def func (x) :
return 1. + np.sin (x) ;
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 6 :
print ('usage: ', sys.argv[0], 'xMin xMax yMin yMax N_evt_min N_evt_max')
exit ()
xMin = float (sys.argv[1])
xMax = float (sys.argv[2])
yMax = float (sys.argv[3])
N_evt_min = int (sys.argv[4])
N_evt_max = int (sys.argv[5])
N_events = []
HOM_integrals = []
HOM_integral_errors = []
CRU_integrals = []
CRU_integral_errors = []
while (N_evt_min < N_evt_max) :
HOM_integral, HOM_integral_unc = integral_HOM (func, xMin, xMax, yMax, N_evt_min)
HOM_integrals.append (HOM_integral)
HOM_integral_errors.append (HOM_integral_unc)
CRU_integral, CRU_integral_unc = integral_CrudeMC (func, xMin, xMax, N_evt_min)
CRU_integrals.append (CRU_integral)
CRU_integral_errors.append (CRU_integral_unc)
N_events.append (N_evt_min)
N_evt_min = N_evt_min * 2
fig, ax = plt.subplots ()
ax.set_title ('integral value and precision', size=14)
ax.set_xlabel ('number of events')
ax.set_ylabel ('integral and its error')
ax.errorbar (N_events, HOM_integrals, xerr = 0.0, yerr = HOM_integral_errors,
color = 'red', label = 'hit or miss')
ax.errorbar (N_events, CRU_integrals, xerr = 0.0, yerr = CRU_integral_errors,
color = 'blue', label = 'crude MC')
ax.legend ()
plt.savefig ('es_8.8.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
8.5.9. Exercise 8.9#
es_8.9.py
#!/usr/bin/python
'''
Si utilizzi il metodo hit-or-miss per stimare l’integrale sotteso
ad una disrtibuzione di probabilita' Gaussiana con *μ=0* e *σ=1
in un generico intervallo *[a,b]*.
* Si calcoli l'integrale contenuto entro gli intervalli *[-kσ,kσ]*
al variare di k da *1* a *5*.
python3 es_8.10.py 100000
'''
import sys
import numpy as np
from integral import integral_HOM
def func (x) :
'''
normal distribution: Gaussian with mean and sigma equal to 1
'''
return 0.3989422804 * np.exp (-0.5 * x * x)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 2 :
print ('usage: ', sys.argv[0], 'N_evt')
exit ()
N_evt = int (sys.argv[1])
for i in range (1, 6) :
integral, integral_unc = integral_HOM (func, float (-i), float (i), func (0), N_evt)
print ('integral at ', i, 'sigma: ', integral,
' +- ', integral_unc)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()