In [1]:
import numpy as np
import matplotlib.pyplot as plt
simulace náhodné proměnné s rozdělením exp + gauss
In [4]:
def sim():
branching_ratio=0.7
r=np.random.random()
if r<branching_ratio:
return np.random.normal(8,1)
else:
return np.random.exponential(3)
In [6]:
nbins=100 #počet binů
ibin=10 #budeme sledovat počty hodnot v 10-tém binu
Nrepeat=100 #počet opakování
Ntot=np.arange(100,5000,100) #celkový počet hodnot v histogramu
value=np.empty(Nrepeat) #pole hodnot v 10-tém binu
fraction=np.empty(np.size(Ntot)) #frakce hodnot připadající na 10-tý bin
std=np.empty(np.size(Ntot)) #standardní odchylky
In [8]:
ik=0
for k in Ntot: #cyklus pro postupně narůstající počet dat Ntot
print('počet dat = ',k)
for j in range(Nrepeat): #Nrepeat opakování pro každé Ntot
data=np.empty(k)
for i in range(k): #generování histogramu
data[i]=sim()
hist,bin_edges=np.histogram(data,bins=nbins)
value[j]=hist[ibin]
fraction[ik]=np.mean(value)/k #jaká frakce hodnot padne do 10-tého binu
std[ik]=np.std(value) #standarní odchylka pro počty v 10-tém binu
ik+=1
fraction_mean=np.mean(fraction) #průměrná frakce připadající na 10-tý bin
počet dat = 100 počet dat = 200 počet dat = 300 počet dat = 400 počet dat = 500 počet dat = 600 počet dat = 700 počet dat = 800 počet dat = 900 počet dat = 1000 počet dat = 1100 počet dat = 1200 počet dat = 1300 počet dat = 1400 počet dat = 1500 počet dat = 1600 počet dat = 1700 počet dat = 1800 počet dat = 1900 počet dat = 2000 počet dat = 2100 počet dat = 2200 počet dat = 2300 počet dat = 2400 počet dat = 2500 počet dat = 2600 počet dat = 2700 počet dat = 2800 počet dat = 2900 počet dat = 3000 počet dat = 3100 počet dat = 3200 počet dat = 3300 počet dat = 3400 počet dat = 3500 počet dat = 3600 počet dat = 3700 počet dat = 3800 počet dat = 3900 počet dat = 4000 počet dat = 4100 počet dat = 4200 počet dat = 4300 počet dat = 4400 počet dat = 4500 počet dat = 4600 počet dat = 4700 počet dat = 4800 počet dat = 4900
In [16]:
plt.scatter(Ntot,std)
plt.plot(Ntot,np.sqrt(Ntot*fraction_mean),c='red')
plt.xlabel('$N_{tot}$',fontsize=14)
plt.ylabel('$\sigma$',fontsize=14)
plt.show()
In [ ]: