Centrální limitní věta

sčítání rovnom. rozdělených dat

  • součet k=5,10,...,25 vzorků z R(0,10) - 1e4 opakování
  • porovnání s N(5 k,sqrt(k*100/12.))

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]:
<matplotlib.axes.AxesSubplot at 0xa4434ac>

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]:
<matplotlib.axes.AxesSubplot at 0x96d1a4c>

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]:
<matplotlib.legend.Legend at 0x983134c>

Maximální hodnota souboru

Soubor o 10,30,100,200,300 prvcích s norm. rozdělením kolem 0


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]:
<matplotlib.axes.AxesSubplot at 0x9e6684c>

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]:
([0.417262945390947,
  0.5441148472334048,
  0.7068838411589191,
  0.7786156480734037,
  0.7654728575262123],
 [0.36461010507330105,
  0.48725063737632324,
  1.0282534723415635,
  1.1808571651136885,
  1.072496167603334])

In [6]:
from numpy import random
random.exponential(2,size=100).max()


Out[6]:
8.5979750780631647

In [7]:
his1[1]


Out[7]:
array([ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,
        0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.0001,
        0.0006,  0.0007,  0.0004,  0.0007,  0.0014,  0.0022,  0.0024,
        0.0039,  0.0056,  0.0081,  0.0095,  0.0122,  0.0155,  0.0177,
        0.0211,  0.0238,  0.0276,  0.0311,  0.0333,  0.0378,  0.0398,
        0.0417,  0.0399,  0.0425,  0.0422,  0.0393,  0.0416,  0.0411,
        0.0362,  0.0343,  0.0333,  0.0322,  0.0325,  0.0267,  0.0241,
        0.0252,  0.0207,  0.0179,  0.0169,  0.0145,  0.0118,  0.0124,
        0.0097,  0.0091,  0.0099,  0.0083,  0.0057,  0.0055,  0.0039,
        0.0043,  0.0032,  0.0027,  0.0026,  0.0025,  0.002 ,  0.0013,
        0.0009,  0.0009,  0.0008,  0.0009,  0.0003,  0.0006,  0.0004,
        0.0005,  0.0004,  0.0003,  0.0001,  0.0001,  0.0001,  0.    ,
        0.0001,  0.0001,  0.0002,  0.0001,  0.    ,  0.    ,  0.    ,
        0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ])