from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'svg')
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()`
%matplotlib inline
import numpy as np
import math
import sys
import random
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
from scipy.stats import gamma
from scipy.stats import poisson
from scipy.stats import binom
from scipy import stats
def MisuraConteggi(T, n_misure, mode='tufo', Rt=1, Rc=5):
# risoluzione temporale del contatore
res = 0.001
# possibilità di combinare sorgenti diverse
if mode == 'tufo':
R = Rt + Rc
elif mode == 'fondo':
R = Rc
else:
print ("Opzione errata! Usa tufo (per misura sul blocco di tufo) o 'fondo' (per misura solo fondo cosmico)")
sys.exit()
print('La rate vera R = ', R)
# lo strumento accetta solo tempi al secondo tra 1 e 600
T = round(T,1)
if T<1 or T>600:
print ("La durata di una singola presa dati deve essere un numero tra 1s e 600s")
sys.exit()
# liste per registrare i conteggi e i tempi di ciascun conteggio
counts_list = []
time_counts_list = []
# assumo l'intevallo pari alla risoluzione
# calcolo il numero di intervalli
n = int(T / res)
# calcolo la probabilita' di conteggio nell'intevallino
p = R * res
if (p > 0.01):
print('Lo strumento non ha una risoluzione temporale per osservare eventi con rate cosi elevata')
sys.exit()
# faccio un loop sul numero di misure
zeri = 0
uni = 0
for mis in range(1, n_misure+1):
# azzero la lista dei tempi ed il contatore di eventi
tmp_time_counts_list = []
tmp_count = 0
# faccio il loop sugli intervallini
for i in range(1, n+1):
# se l'evento accade, registro il tempo come:
rrr = random.random()
if (rrr < p):
# opportuno indice che rappresenta l'intervallino
# in cui accade un evento
hit_time = round(i*res, 4) # misuro il tempo come
# i * Delta t (in questo caso la risoluzione)
tmp_time_counts_list.append(hit_time) # metto il tempo dell'evento in
# esco dal loop sugli intervallini
# conto gli eventi, e copio i risutati in una lista contenente
# le varie misure
counts_list.append(len(tmp_time_counts_list))
if (len(tmp_time_counts_list) == 0):
zeri += 1
if (len(tmp_time_counts_list) == 1):
uni += 1
# copio tutti i tempi degli eventi (un array) in una lista contentente
# le varie misure
time_counts_list.append(np.asarray(tmp_time_counts_list))
### esco dal loop sulle misure
################################################
# converto tutto in numpy-arrays
time_counts_np = np.asarray(time_counts_list)
counts_np = np.asarray(counts_list)
myzeri= zeri/n_misure
myuni=(zeri+uni)/n_misure
myzeri_err= math.sqrt(myzeri*(1-myzeri)/n_misure)
myuni_err= math.sqrt(myuni*(1-myuni)/n_misure)
if (myzeri_err < 0.05):
myzeri_err = 0.05
if (myuni_err < 0.05):
myuni_err = 0.05
print ("{0:d} misure ripetute di conteggio in un intervallo temporale di {1:.1f} s =".format(n_misure,T))
print (counts_np)
print ("Frazione di conteggi di zero:", myzeri,"+-", myzeri_err)
print ("Frazione di conteggi di zero ed uno:", myuni, "+-", myuni_err)
#print (R)
return counts_np, time_counts_np, myzeri, myzeri_err, myuni,myuni_err
def exp(t, R):
return R * np.exp(-R*t)
def exp_2(t, R):
return t*R*R*np.exp(-R*t)
def tempoattesa(mytimes, k):
t_list = []
if(k==1):
for i in range (0, mytimes.size):
if(mytimes[i].size > 1):
for j in range (mytimes[i].size-1, 0, -1):
diff = mytimes[i][j] - mytimes[i][j-1]
t_list.append(diff)
t = np.asarray(t_list)
if(k==2):
for i in range (0, mytimes.size):
if(mytimes[i].size > 2):
for j in range (mytimes[i].size-1, 1, -1):
diff = mytimes[i][j] - mytimes[i][j-2]
t_list.append(diff)
t = np.asarray(t_list)
return t
random.seed(6516654)
dt_1 = 1
misure = 50
T_1 = dt_1*misure
mycounts_1, mytimes_1, myzeri_1, myzeri_err_1, myuni_1, myuni_err_1 = MisuraConteggi(dt_1, misure, mode= 'fondo', Rc=1.21, Rt=0.60)
x_1 = mycounts_1.sum()
lambda_1 = (x_1+1)/misure
sigma_lambda_1 = math.sqrt(lambda_1)/math.sqrt(misure)
R_1 = lambda_1/dt_1
sigma_R_1 = sigma_lambda_1/dt_1
print ("lambda =", lambda_1)
print ("sigma_lambda =", sigma_lambda_1)
print ("R =", R_1)
print ("sigma_R =", sigma_R_1)
50 misure ripetute di conteggio in un intervallo temporale di 1.0 s =
[0 4 2 1 1 2 1 1 1 0 5 1 2 3 2 1 0 1 1 2 1 3 1 1 0 0 0 2 0 2 2 3 0 2 1 1 3
0 5 1 0 2 1 1 1 2 1 0 1 3]
Frazione di conteggi di zero: 0.22 +- 0.05858327406350724
Frazione di conteggi di zero ed uno: 0.62 +- 0.06864400920692205
lambda = 1.44
sigma_lambda = 0.1697056274847714
R = 1.44
sigma_R = 0.1697056274847714
/shared-libs/python3.7/py/lib/python3.7/site-packages/numpy/core/_asarray.py:83: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
return array(a, dtype, copy=False, order=order)
dt_2 = 2
misure = 50
T_2 = dt_2*misure
mycounts_2, mytimes_2, myzeri_2, myzeri_err_2, myuni_2, myuni_err_2 = MisuraConteggi(dt_2, misure, mode= 'fondo', Rc=1.21, Rt=0.60)
x_2 = mycounts_2.sum()
lambda_2 = (x_2+1)/misure
sigma_lambda_2 = math.sqrt(lambda_2)/math.sqrt(misure)
R_2 = lambda_2/dt_2
sigma_R_2 = sigma_lambda_2/dt_2
print ("lambda =", lambda_2)
print ("sigma_lambda =", sigma_lambda_2)
print ("R =", R_2)
print ("sigma_R =", sigma_R_2)
50 misure ripetute di conteggio in un intervallo temporale di 2.0 s =
[3 4 0 3 3 2 2 2 1 2 0 2 2 1 4 1 3 5 3 1 2 7 3 0 5 1 0 1 4 6 2 3 2 1 4 5 3
1 2 3 5 5 0 3 2 5 2 4 1 1]
Frazione di conteggi di zero: 0.1 +- 0.05
Frazione di conteggi di zero ed uno: 0.3 +- 0.0648074069840786
lambda = 2.56
sigma_lambda = 0.2262741699796952
R = 1.28
sigma_R = 0.1131370849898476
dt_3 = 3
misure = 50
T_3 = dt_3*misure
mycounts_3, mytimes_3, myzeri_3, myzeri_err_3, myuni_3, myuni_err_3 = MisuraConteggi(dt_3, misure, mode= 'fondo', Rc=1.21, Rt=0.60)
x_3 = mycounts_3.sum()
lambda_3 = (x_3+1)/misure
sigma_lambda_3 = math.sqrt(lambda_3)/math.sqrt(misure)
R_3 = lambda_3/dt_3
sigma_R_3 = sigma_lambda_3/dt_3
print ("lambda =", lambda_3)
print ("sigma_lambda =", sigma_lambda_3)
print ("R =", R_3)
print ("sigma_R =", sigma_R_3)
50 misure ripetute di conteggio in un intervallo temporale di 3.0 s =
[3 2 3 2 6 3 5 1 2 3 7 4 2 6 2 8 2 1 2 2 2 2 6 4 4 4 9 4 6 1 4 3 7 4 6 3 3
2 4 5 4 1 6 3 4 4 3 5 4 3]
Frazione di conteggi di zero: 0.0 +- 0.05
Frazione di conteggi di zero ed uno: 0.08 +- 0.05
lambda = 3.74
sigma_lambda = 0.27349588662354685
R = 1.2466666666666668
sigma_R = 0.09116529554118229
dt_4 = 4
misure = 50
T_4 = dt_4*misure
mycounts_4, mytimes_4, myzeri_4, myzeri_err_4, myuni_4, myuni_err_4 = MisuraConteggi(dt_4, misure, mode= 'fondo', Rc=1.21, Rt=0.60)
x_4 = mycounts_4.sum()
lambda_4 = (x_4+1)/misure
sigma_lambda_4 = math.sqrt(lambda_4)/math.sqrt(misure)
R_4 = lambda_4/dt_4
sigma_R_4 = sigma_lambda_4/dt_4
print ("lambda =", lambda_4)
print ("sigma_lambda =", sigma_lambda_4)
print ("R =", R_4)
print ("sigma_R =", sigma_R_4)
50 misure ripetute di conteggio in un intervallo temporale di 4.0 s =
[5 1 2 6 5 8 6 2 6 4 7 1 3 6 2 8 4 6 4 4 3 5 6 2 5 6 5 3 4 5 5 4 3 6 2 5 5
5 4 2 4 4 3 5 7 4 6 7 8 6]
Frazione di conteggi di zero: 0.0 +- 0.05
Frazione di conteggi di zero ed uno: 0.04 +- 0.05
lambda = 4.6
sigma_lambda = 0.30331501776206204
R = 1.15
sigma_R = 0.07582875444051551
dt_10 = 10
misure = 50
T_10 = dt_10*misure
mycounts_10, mytimes_10, myzeri_10, myzeri_err_10, myuni_10, myuni_err_10 = MisuraConteggi(dt_10, misure, mode= 'fondo', Rc=1.21, Rt=0.60)
x_10 = mycounts_10.sum()
lambda_10 = (x_10+1)/misure
sigma_lambda_10 = math.sqrt(lambda_10)/math.sqrt(misure)
R_10 = lambda_10/dt_10
sigma_R_10 = sigma_lambda_10/dt_10
print ("lambda =", lambda_10)
print ("sigma_lambda =", sigma_lambda_10)
print ("R =", R_10)
print ("sigma_R =", sigma_R_10)
50 misure ripetute di conteggio in un intervallo temporale di 10.0 s =
[12 17 15 8 14 10 16 10 14 8 15 10 7 8 17 11 16 13 14 16 9 10 17 13
10 17 11 21 13 9 12 9 11 10 9 10 11 17 13 14 13 16 10 15 14 11 11 11
12 5]
Frazione di conteggi di zero: 0.0 +- 0.05
Frazione di conteggi di zero ed uno: 0.0 +- 0.05
lambda = 12.32
sigma_lambda = 0.4963869458396343
R = 1.232
sigma_R = 0.04963869458396343
dt = 599
misure = 1
T = dt*misure
mycounts, mytimes, myzeri, myzeri_err, myuni, myuni_err = MisuraConteggi(dt, misure, mode= 'fondo', Rc=1.21, Rt=0.60)
x = mycounts.sum()
lambda_ = (x+1)/misure
sigma_lambda = math.sqrt(lambda_)/math.sqrt(misure)
R = lambda_/dt
sigma_R = sigma_lambda/dt
print ("lambda =", lambda_)
print ("sigma_lambda =", sigma_lambda)
print ("R =", R)
print ("sigma_R =", sigma_R)
1 misure ripetute di conteggio in un intervallo temporale di 599.0 s =
[751]
Frazione di conteggi di zero: 0.0 +- 0.05
Frazione di conteggi di zero ed uno: 0.0 +- 0.05
lambda = 752.0
sigma_lambda = 27.422618401604176
R = 1.2554257095158599
sigma_R = 0.04578066511119228
x_tot = x_1 + x_2 + x_3 + x_4 + x_10 + x
T_tot = T_1 + T_2 + T_3 + T_4 + T_10 + T
R_best = (x_tot+1)/T_tot
sigma_R_best = math.sqrt(x_tot+1)/T_tot
R_best
sigma_R_best
Istogrammi
bins = np.arange(-0.5, mycounts_1.max()+0.5)
prova1 = plt.hist(mycounts_1, bins=bins, density=True)
plt.ylabel('Frequenza')
plt.xlabel('Numero di conteggi')
plt.title('50 misure da 1s')
bins = np.arange(-0.5, mycounts_2.max()+0.5)
plt.hist(mycounts_2, bins=bins, density=True)
plt.ylabel('Frequenza')
plt.xlabel('Numero di conteggi')
plt.title('50 misure da 2s')
bins = np.arange(-0.5, mycounts_3.max()+0.5)
plt.hist(mycounts_3, bins=bins, density=True)
plt.ylabel('Frequenza')
plt.xlabel('Numero di conteggi')
plt.title('50 misure da 3s')
bins = np.arange(-0.5, mycounts_4.max()+0.5)
plt.hist(mycounts_4, bins=bins, density=True)
plt.ylabel('Frequenza')
plt.xlabel('Numero di conteggi')
plt.title('50 misure da 4s')
bins = np.arange(-0.5, mycounts_10.max()+0.5)
plt.hist(mycounts_10, bins=bins,density=True)
plt.ylabel('Frequenza')
plt.xlabel('Numero di conteggi')
plt.title('50 misure da 10s')
mean_1, var_1 = gamma.stats(a=(x_1+1), moments='mv')
print('media t1 = ', mean_1/T_1)
print('deviazione standard = ', np.sqrt(var_1)/T_1)
mean_4, var_4 = gamma.stats(a=(x_4+1), moments='mv')
print('media t4 = ', mean_4/T_4)
print('deviazione standard = ', np.sqrt(var_4)/T_4)
mean, var = gamma.stats(a=(x+1), moments='mv')
print('media t599 = ', mean/T)
print('deviazione standard = ', np.sqrt(var)/T)
mean_tot, var_tot = gamma.stats(a=(x_tot+1), moments='mv')
print('media combinata= ', mean_tot/T_tot)
print('deviazione standard = ', np.sqrt(var_tot)/T_tot)
estremi =np.array([mean_1/T_1-(3*np.sqrt(var_1)/T_1), mean_1/T_1+(3*np.sqrt(var_1)/T_1),
mean_4/T_4-(3*np.sqrt(var_4)/T_4), mean_4/T_4+(3*np.sqrt(var_4)/T_4),
mean/T-3*np.sqrt(var)/T, mean/T+3*np.sqrt(var)/T,
mean_tot/T_tot-3*np.sqrt(var_tot)/T_tot, mean_tot/T_tot+3*np.sqrt(var_tot)/T_tot])
r = np.linspace(estremi.min(), estremi.max(), 1000)
graph_1= gamma.pdf(r*T_1, a=(x_1+1))*T_1
graph_4= gamma.pdf(r*T_4, a=(x_4+1))*T_4
graph= gamma.pdf(r*T, a=(x+1))*T
graph_tot= gamma.pdf(r*(T_tot), a=(x_tot+1))*T_tot
mode_1 = r[np.argmax(graph_1)]
mode_4 = r[np.argmax(graph_4)]
mode = r[np.argmax(graph)]
mode_tot= r[np.argmax(graph_tot)]
print('moda t1 =', mode_1)
print('moda t4 =', mode_4)
print('moda t599 =', mode)
print('moda combinata =', mode_tot)
figure(figsize=(10, 5), dpi=80)
plt.vlines(1.21, 0, 15,linestyles="--", label='Valore vero')
plt.plot(r, gamma.pdf(r*T_1, a=(x_1+1))*T_1, label = '$R_1$')
plt.plot(r, gamma.pdf(r*T_4, a=(x_4+1))*T_4, label = '$R_4$')
plt.plot(r, gamma.pdf(r*T, a=(x+1))*T, label = '$R_{599}$')
plt.plot(r, gamma.pdf(r*(T_tot), a=(x_tot+1))*T_tot, label = '$R_{best}$')
plt.xlabel('Tasso di conteggio $R$', fontsize=16)
plt.ylabel('$pdf(R)$', fontsize=16)
plt.title('Posterior pdf del tasso di conteggio', fontsize=16)
plt.legend()
media t1 = 1.44
deviazione standard = 0.1697056274847714
media t4 = 1.15
deviazione standard = 0.07582875444051551
media t599 = 1.2554257095158599
deviazione standard = 0.04578066511119228
media combinata= 1.2382739212007505
deviazione standard = 0.02782813691217885
moda t1 = 1.4198870325298214
moda t4 = 1.145509615190327
moda t599 = 1.2534108467283305
moda combinata = 1.2379963850800442
Punto 5
l = R_best*dt_1
p=poisson(mu=l)
#n = l*2+3
n = mycounts_1.max()
bins = np.arange(-0.5+mycounts_1.min(), mycounts_1.max()+0.5)
plt.hist(mycounts_1, bins=bins,density=True, label='Dati osservati')
x = np.arange(mycounts_1.min(), n)
# grafico della pmf
plt.plot(x, p.pmf(x), 'bo', ms=4, color='b', label='Distribuzione Poisson con $R_{best}$')
plt.vlines(x, 0, p.pmf(x), colors='b', lw=2, alpha=0.5)
plt.title('Conteggi per 50 misure da 1s', fontsize=16)
plt.xlabel("Numero di conteggi", fontsize=16)
plt.ylabel("Frequenza", fontsize=16)
plt.legend()
t = plt.xticks(x)
#plt.ylim(0, np.max(p.pmf(x))*2)
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:11: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "bo" (-> color='b'). The keyword argument will take precedence.
# This is added back by InteractiveShellApp.init_path()
l = R_best*dt_2
p=poisson(mu=l)
#n = l*2+3
n = mycounts_2.max()
bins = np.arange(-0.5+mycounts_2.min(), mycounts_2.max()+0.5)
plt.hist(mycounts_2, bins=bins,density=True, label='Dati osservati')
x = np.arange(mycounts_2.min(), n)
# grafico della pmf
plt.plot(x, p.pmf(x), 'bo', ms=4, label='Distribuzione Poisson con $R_{best}$')
plt.vlines(x, 0, p.pmf(x), colors='b', lw=2, alpha=0.5)
plt.title('Conteggi per 50 misure da 2s', fontsize=16)
plt.xlabel("Numero di conteggi", fontsize=16)
plt.ylabel("Frequenza", fontsize=16)
plt.legend()
t = plt.xticks(x)
#plt.ylim(0, np.max(p.pmf(x))*2)
l = R_best*dt_3
p=poisson(mu=l)
#n = l*2+3
n = mycounts_3.max()
bins = np.arange(-0.5+mycounts_3.min(), mycounts_3.max()+0.5)
plt.hist(mycounts_3, bins=bins,density=True, label='Dati osservati')
x = np.arange(mycounts_3.min(), n)
# grafico della pmf
plt.plot(x, p.pmf(x), 'bo', ms=4, label='Distribuzione Poisson con $R_{best}$')
plt.vlines(x, 0, p.pmf(x), colors='b', lw=2, alpha=0.5)
plt.title('Conteggi per 50 misure da 3s', fontsize=16)
plt.xlabel("Numero di conteggi", fontsize=16)
plt.ylabel("Frequenza", fontsize=16)
plt.legend()
t = plt.xticks(x)
#plt.ylim(0, np.max(p.pmf(x))*2)
l = R_best*dt_4
p=poisson(mu=l)
#n = l*2+3
n = mycounts_4.max()
bins = np.arange(-0.5+mycounts_4.min(), mycounts_4.max()+0.5)
plt.hist(mycounts_4, bins=bins,density=True, label='Dati osservati')
x = np.arange(mycounts_4.min(), n)
# grafico della pmf
plt.plot(x, p.pmf(x), 'bo', ms=4, label='Distribuzione Poisson con $R_{best}$')
plt.vlines(x, 0, p.pmf(x), colors=