In [8]:
from math import sqrt,pi
from numpy import random,histogram,r_,array,exp,polyfit
nexp=10000
sim=array([random.uniform(0,10,30).reshape(5,6).sum(0).cumsum() for i in range(nexp)]).transpose()
bin=r_[0:200:1.]
import matplotlib.pyplot as rp
his1=[histogram(sim[0],bin)[0]/float(nexp)]
a=rp.plot(bin[:-1],his1[0])[0].axes
i=0
sig=sqrt((i+1)*500./12.)
a.plot(bin[:-1],1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2),'k:')
for i in range(1,len(sim)):
s=sim[i]
his1.append(histogram(s,bin)[0]/float(nexp))
sig=sqrt((i+1)*500./12.)
a.plot(bin[:-1],his1[-1])
a.plot(bin[:-1],1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2),'k:')
sum(his1[-3])
a.legend(r_[5:31:5])
a
Out[8]:
exponenciální rozdělení - E(5) o stejných počtech, srovnání s N(5k,25k)
In [2]:
sim=array([random.exponential(5,30).reshape(5,6).sum(0).cumsum() for i in range(nexp)]).transpose()
bin=r_[0:200:1.]
import matplotlib.pyplot as rp
his1=[histogram(sim[0],bin)[0]/float(nexp)]
a=rp.plot(bin[:-1],his1[0])[0].axes
i=0
sig=sqrt((i+1)*125.)
a.plot(bin[:-1],1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2),'k:')
for i in range(1,len(sim)):
s=sim[i]
his1.append(histogram(s,bin)[0]/float(nexp))
sig=sqrt((i+1)*125.)
a.plot(bin[:-1],his1[-1])
a.plot(bin[:-1],1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2),'k:')
sum(his1[-3])
a.legend(r_[5:30:5])
a
Out[2]:
In [3]:
for i in range(1,len(sim)):
sig=sqrt((i+1)*125.)
plot(bin[:-1],his1[i]-1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2))
legend(r_[5:30:5])
Out[3]:
In [4]:
psim2=[]
bs=[10,30,100,200,300]
for i in range(nexp):
dat=random.normal(0,1.,300)
psim2.append([max(dat[:k]) for k in bs])
import matplotlib.pyplot as rp
sim2=array(psim2).transpose()
bin=r_[:5:0.05]
his1=[histogram(sim2[0],bin)[0]/float(nexp)]
b=rp.plot(bin[:-1],his1[0])[0].axes
i=0
sig=sqrt((i+1)*500./12.)
#b.plot(bin[:-1],1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2),'k:')
for i in range(1,len(sim2)):
s=sim2[i]
his1.append(histogram(s,bin)[0]/float(nexp))
sig=sqrt((i+1)*500./12.)
b.plot(bin[:-1],his1[-1])
#b.plot(bin[:-1],1/sqrt(2*pi)/sig*exp(-(bin[:-1]-(i+1)*25)**2/2/sig**2),'k:')
b.legend(bs)
b.set_xlabel("max.hodnota opak. norm. rozdeleni")
b
Out[4]:
použití vyšších momentů (šikmost, špičatost) pro testy normálnosti (viz)
In [5]:
from scipy import stats
ta=[stats.skew(h) for h in sim2]
tb=[stats.kurtosis(h) for h in sim2]
ta,tb
Out[5]:
In [6]:
from numpy import random
random.exponential(2,size=100).max()
Out[6]:
In [7]:
his1[1]
Out[7]: