7.6. Solutions: Lezione 7#
7.6.1. Exercise 7.1#
es_7.1.py
#!/usr/bin/python
'''
Si generi un campione di numeri pseudo-casuali
distribuiti secondo una distribuzione di densità esponenziale
con tempo caratteristico t<sub>0</sub> di 5 secondi
e si visualizzi la distribuzione del campione ottenuto
in un istogramma utilizzando il metodo della funzione inversa.
Si scrivano tutte le funzioni deputate alla generazione di numeri casuali
in una libreria, implementata in file separati rispetto al programma principale.
python3 es_7.1.py 1000 3.1 15
'''
import sys
import matplotlib.pyplot as plt
import numpy as np
from math import floor
from myrand import generate_exp
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 3 :
print ('usage: ', sys.argv[0], 'N tau [seed=10]')
exit ()
tau = float (sys.argv[2])
if tau <= 0 :
print ('The tau parameter has to be positive')
exit ()
N = int (sys.argv[1])
seed = 10
if len(sys.argv) >= 4 :
seed = float (sys.argv[3])
randlist = generate_exp (tau, N, seed)
# plotting of the generated list of numbers in a histogram
nBins = floor (len (randlist) / 100.) # number of bins of the hitogram
bin_edges = np.linspace (0., 4. * tau, nBins + 1) # edges o the histogram bins
fig, ax = plt.subplots ()
ax.set_title ('Exponential pdf with tau ' + str (tau), size=14)
ax.set_xlabel ('random value')
ax.set_ylabel ('events in bin')
ax.hist (randlist, # list of numbers
bins = bin_edges,
color = 'orange',
# normed = True,
)
plt.savefig ('es_7.1.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
7.6.2. Exercise 7.2#
es_7.2.py
#!/usr/bin/python
'''
Si utilizzi il risultato del primo esercizio per simulare uno pseudo-esperimento di conteggio
con caratteristiche di Poisson:
* si scelga un tempo caratteristico t<sub>0</sub> di un processo di decadimento radioattivo;
* si scelta un tempo di misura t<sub>M</sub> entro cui fare conteggi;
* in un ciclo, si simulino N pseudo-esperimenti di conteggio, in cui,
per ciascuno di essi, si simuli una sequenza di eventi casuali
con intertempo caratteristico dei fenomeni di Poisson,
fino a che il tempo totale trascorso non sia maggiore del tempo di misura,
contando il numero di eventi generati che cascano nell'intervallo;
* si riempia un istogramma con i conteggi simulati per ogni esperimento
python3 es_7.2.py 1000 2.2 6.6
'''
import sys
import matplotlib.pyplot as plt
import numpy as np
from math import floor, ceil
from myrand import rand_exp
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 4 :
print ('usage: ', sys.argv[0], 'N tau_gen tau_meas [seed=10]')
exit ()
tau_gen = float (sys.argv[2])
if tau_gen <= 0 :
print ('The tau_gen parameter has to be positive')
exit ()
tau_meas = float (sys.argv[3])
if tau_meas < tau_gen :
print ('The tau_meas parameter has to be larger than tau_gen')
exit ()
N = int (sys.argv[1])
seed = 10
if len(sys.argv) >= 5 :
seed = float (sys.argv[4])
randlist = []
for i in range (N) :
total_time = rand_exp (tau_gen)
events = 0
while (total_time < tau_meas) :
events = events + 1
total_time = total_time + rand_exp (tau_gen)
randlist.append (events)
my_stats = stats (randlist)
# plotting of the generated list of numbers in a histogram
nBins = floor (len (randlist) / 100.) # number of bins of the hitogram
xMin = max (0., ceil (my_stats.mean () - 3 * my_stats.sigma ()))
xMax = ceil (my_stats.mean () + 3 * my_stats.sigma ())
bin_edges = np.linspace (xMin, xMax, int (xMax - xMin) + 1) # edges o the histogram bins
fig, ax = plt.subplots ()
ax.set_title ('Counts distribution', size=14)
ax.set_xlabel ('random value')
ax.set_ylabel ('events in bin')
ax.hist (randlist, # list of numbers
bins = bin_edges,
color = 'orange',
# normed = True,
)
plt.savefig ('es_7.2.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
7.6.3. Exercise 7.3#
es_7.3.cpp
#!/usr/bin/python
'''
Si utilizzi il risultato del primo esercizio per simulare uno pseudo-esperimento di conteggio
con caratteristiche di Poisson:
* si scelga un tempo caratteristico t<sub>0</sub> di un processo di decadimento radioattivo;
* si scelta un tempo di misura t<sub>M</sub> entro cui fare conteggi;
* in un ciclo, si simulino N pseudo-esperimenti di conteggio, in cui,
per ciascuno di essi, si simuli una sequenza di eventi casuali
con intertempo caratteristico dei fenomeni di Poisson,
fino a che il tempo totale trascorso non sia maggiore del tempo di misura,
contando il numero di eventi generati che cascano nell'intervallo;
* si riempia un istogramma con i conteggi simulati per ogni esperimento
python3 es_7.3.py 1000 4.2
'''
import sys
import matplotlib.pyplot as plt
import numpy as np
from math import floor, ceil
from myrand import generate_poisson
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 3 :
print ('usage: ', sys.argv[0], 'N mean [seed=10]')
exit ()
mean = float (sys.argv[2])
if mean <= 0 :
print ('The mean parameter has to be positive')
exit ()
N = int (sys.argv[1])
seed = 10
if len(sys.argv) >= 4 :
seed = float (sys.argv[3])
randlist = generate_poisson (mean, N)
my_stats = stats (randlist)
print ('mean : ', mean, my_stats.mean ())
print ('variance : ', mean, my_stats.variance ())
# plotting of the generated list of numbers in a histogram
nBins = floor (len (randlist) / 100.) # number of bins of the hitogram
xMin = max (0., ceil (my_stats.mean () - 3 * my_stats.sigma ()))
xMax = ceil (my_stats.mean () + 3 * my_stats.sigma ())
bin_edges = np.linspace (xMin, xMax, int (xMax - xMin) + 1) # edges o the histogram bins
fig, ax = plt.subplots ()
ax.set_title ('Counts distribution', size=14)
ax.set_xlabel ('random value')
ax.set_ylabel ('events in bin')
ax.hist (randlist, # list of numbers
bins = bin_edges,
color = 'orange',
# normed = True,
)
plt.savefig ('es_7.3.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()
stats.py
#!/usr/bin/python
from math import sqrt, pow
class stats :
'''calculator for statistics of a list of numbers'''
summ = 0.
sumSq = 0.
N = 0
sample = []
def __init__ (self, sample):
'''
reads as input the collection of events,
which needs to be a list of numbers
'''
self.sample = sample
self.summ = sum (self.sample)
self.sumSq = sum ([x*x for x in self.sample])
self.N = len (self.sample)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def mean (self) :
'''
calculates the mean of the sample present in the object
'''
return self.summ / self.N
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def variance (self, bessel = True) :
'''
calculates the variance of the sample present in the object
'''
var = self.sumSq / self.N - self.mean () * self.mean ()
if bessel : var = self.N * var / (self.N - 1)
return var
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def sigma (self, bessel = True) :
'''
calculates the sigma of the sample present in the object
'''
return sqrt (self.variance (bessel))
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def sigma_mean (self, bessel = True) :
'''
calculates the sigma of the sample present in the object
'''
return sqrt (self.variance (bessel) / self.N)
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def skewness (self) :
'''
calculate the skewness of the sample present in the object
'''
mean = self.mean ()
asymm = 0.
for x in self.sample:
asymm = asymm + pow (x - mean, 3)
asymm = asymm / (self.N * pow (self.sigma (), 3))
return asymm
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def kurtosis (self) :
'''
calculate the kurtosis of the sample present in the object
'''
mean = self.mean ()
kurt = 0.
for x in self.sample:
kurt = kurt + pow (x - mean, 4)
kurt = kurt / (self.N * pow (self.variance (), 2)) - 3
return kurt
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
def append (self, x):
'''
add an element to the sample
'''
self.sample.append (x)
self.summ = self.summ + x
self.sumSq = self.sumSq + x * x
self.N = self.N + 1
7.6.4. Exercise 7.4#
es_7.4.py
#!/usr/bin/python
'''
* Si utilizzi il risultato dell'esercizio precedente
per calcolare le statistiche di una distribuzione di Poisson
al variare della media, fra 1 e 250 (come conviene campionare l'intervallo?).
* Si disegni l'andamento ottenuto della skewness e della curtosi in funzione della media
python3 es_7.4.py 10000 1 500
'''
import sys
import matplotlib.pyplot as plt
import numpy as np
from math import floor, ceil
from myrand import generate_poisson
from stats import stats
def main () :
'''
Funzione che implementa il programma principale
'''
if len(sys.argv) < 3 :
print ('usage: ', sys.argv[0], 'N mean_min mean_max [seed=10]')
exit ()
mean_min = float (sys.argv[2])
if mean_min <= 0 :
print ('The mean parameter has to be positive')
exit ()
mean_max = float (sys.argv[3])
if mean_max < mean_min :
print ('The mean_max parameter has to be larger than the mean_min one')
exit ()
N = int (sys.argv[1])
seed = 10
if len(sys.argv) >= 5 :
seed = float (sys.argv[4])
means = []
skewn = []
kurts = []
while (mean_min < mean_max) :
randlist = generate_poisson (mean_min, N)
my_stats = stats (randlist)
means.append (my_stats.mean ())
skewn.append (my_stats.skewness ())
kurts.append (my_stats.kurtosis ())
mean_min *= 2
fig, (ax1, ax2) = plt.subplots (2, sharex=True)
fig.suptitle ('Poisson distribution asymptotic behaviour')
ax1.plot (means, skewn)
ax1.set_ylabel ('skewness')
ax2.plot (means, kurts)
ax2.set_ylabel ('kurtosis')
ax2.set_xlabel ('mean')
plt.savefig ('es_7.4.png')
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
if __name__ == "__main__":
main ()