Ilustrace binomickeho rozdeleni

navazuje na Intervaly spolehlivosti

In [2]:
%pylab inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = [10, 5]
from scipy import stats
n=100
p=0.1
i=30
ptrue=r_[.05:1:0.05]
cont=0.01
xlo=r_[[stats.binom(n,p).ppf(cont) for p in ptrue]]
xhi=r_[[stats.binom(n,p).ppf(1-cont) for p in ptrue]]
plot(xlo,ptrue,':')
plot(xhi,ptrue,':')
cont=0.05
xlo=r_[[stats.binom(n,p).ppf(cont) for p in ptrue]]
xhi=r_[[stats.binom(n,p).ppf(1-cont) for p in ptrue]]
plot(xlo,ptrue)
plot(xhi,ptrue)
plot([0,100],[0,1.])
axvline(20)
xlabel("merena hodnota (100 pokusu)")
ylabel("skutecna hodnota")
legend(["low 1%","high 1%","low 5%","high 5%"])

Populating the interactive namespace from numpy and matplotlib
Out[2]:
<matplotlib.legend.Legend at 0x1ff8d21c128>

(relativní) pološířka binomického rozdělení klesá s počtem pokusů (a blíží se - pro pravděp. blízkou středu 0.5 - normálnímu rozdělení)

In [3]:
plot(r_[:80]/80.,stats.binom(80,0.3).pmf(r_[:80])*4.)
plot(r_[:40]/40.,stats.binom(40,0.3).pmf(r_[:40])*2.)
plot(r_[:20]/20.,stats.binom(20,0.3).pmf(r_[:20]))
stats.binom(80,0.3).pmf(r_[0:40:4])
legend(["%i pokusu"%i for i in [80,40,20]])
Out[3]:
<matplotlib.legend.Legend at 0x1ff8d33ce48>

výpočet mezí pro skutečnou hodnotu podle naměřené dosáhneme "inverzí" funkcí odpovídající modré a zelené mezi pásu spolehlivosti v prvním grafu (konstruujeme horizontálně, odečítáme vertikálně)

In [4]:
#jednoduchá lineární interpolace
interl = lambda a,b,z:b[0]+(b[1]-b[0])*(z-a[0])/(a[1]-a[0])

def truelims(z=20):
    i,j=(xlo<z).sum(),(xhi<z).sum()
    return interl(xhi[j-1:j+1],ptrue[j-1:j+1],z),interl(xlo[i-1:i+1],ptrue[i-1:i+1],z)

truelims(20),truelims(40)
Out[4]:
((0.1416666666666667, 0.27), (0.32, 0.48))
In [10]:
ptrue=r_[.05:1:0.02]
fmap=r_[[stats.binom(80,ptrue[i]).pmf(r_[:80]) for i in range(len(ptrue))]]
imshow(fmap,origin="lower",extent=[0,80,0,1],aspect="auto")
Out[10]:
<matplotlib.image.AxesImage at 0x1ff8d6644a8>
In [11]:
plot(fmap.sum(0)[:70])
plot(fmap.sum(1)[:45])
fmap.sum(0)[:8]
legend(["suma sloupce","suma radku"],loc=4)
Out[11]:
<matplotlib.legend.Legend at 0x1ff8e8c2e80>

pro malé hodnoty je pravděpodob. obsah nízký

Poissonovské rozdělení

In [13]:
#rucni vypocet...
from scipy.special import factorial
poiss=lambda p,n:pow(p,n)/factorial(n)*exp(-p)
mod1=poiss(5.,r_[1:20])
plot(r_[1:20],mod1)
sum(mod1)
Out[13]:
0.9932617077873325
In [14]:
from scipy import optimize
uplim=lambda n,prob,pinit:optimize.fsolve(lambda p:sum(poiss(p,[r_[:n]]))-prob,pinit)[0]
def get_poiss_lims(prob=0.05,nmax=15,pinit=4.):
    rep=[pinit]
    for i in range(1,15):
        rep.append(uplim(i,prob,rep[-1]))
    return rep[1:]
#sum(poiss(4.,r_[:5]))-0.05
In [15]:
nmax=15
l1=get_poiss_lims(0.9,nmax,pinit=1.)
l2=get_poiss_lims(0.95,nmax,pinit=1.)
l3=get_poiss_lims(0.98,nmax,pinit=1.)
u1=get_poiss_lims(0.1,nmax,pinit=4.)
u2=get_poiss_lims(0.05,nmax,pinit=4.)
u3=get_poiss_lims(0.02,nmax,pinit=4.)
In [16]:
x=r_[1:nmax]
[plot(x,y) for y in [l1,l2,l3,u1,u2,u3]]
plot(x,x)
slist=[10,5,2]
legend(["l.l. %i%%"%a for a in slist]+["u.l. %i%%"%a for a in slist],loc=2)
xlabel("mereno")
ylabel("stred. hodnota")
axvline(7)
Out[16]:
<matplotlib.lines.Line2D at 0x1ff8ea044e0>

jednoduchá konstrukce při záměně os

použijeme kvantily Poissonova rozdělení

In [17]:
ktrue=r_[3:25:2]
ktrue2=r_[10:35:2]
for cont in [0.1,.05,0.02]:
    xlo=[stats.poisson(k).ppf(cont) for k in ktrue2]
    xhi=[stats.poisson(k).ppf(1-cont) for k in ktrue]
    fun=np.polyfit(xlo,ktrue2,6)
    fun2=np.polyfit(xhi,ktrue,6)
    plot(xlo,ktrue2,'d')
    plot(xlo,polyval(fun,xlo),'b')
    plot(xhi,polyval(fun2,xhi),'b')
    ylim(0,40)
In [18]:
stats.poisson(15).ppf(0.95)
Out[18]:
22.0

Horni limita

pravděpodobnost v Bayesovském přístupu vychází z poměru (integrálních) pravděpodobností

$$1-\beta=\int_{-\infty}^{\theta_{up}} p(\theta|x) d\theta = \frac{\int_{-\infty}^{\theta_{up}} L(x|\theta) \pi(\theta) d\theta}{\int_{-\infty}^{\infty} L(x|\theta) \pi(\theta) d\theta}$$

kde $\pi(\theta)$ je zvolená a-priorní pravděpodobnost pro danou hodnotu parametru

In [19]:
from scipy import integrate,optimize
cpos=r_[-3:2:0.1]
clasc=stats.norm(cpos,1).isf(0.05)

V nejjednodušším případě $\pi(\theta)$ je konstantní.

Pokud klademe podmínku nezápornosti, pomocí inverse survival function můžeme dopočítat polohu horní limity, kde integrál nad ní je dané procento integrálu přes celou povolenou (nezápornou) oblast.

$L(x|\theta)$ bereme opět normální rozdělení s jednotkovou disperzí.

In [20]:
# cpos je stredni hodnota 
bayes=stats.norm(cpos,1).isf(stats.norm(cpos,1).sf(0)*0.05)
# hodnoty v prom. bayes budou odpovidat ne prsti 0.05, ale 0.05*stats.norm(cpos,1).sf(0)
plot(cpos,clasc)
plot(cpos,bayes)
fill_between(cpos,bayes,4,alpha=0.1)
axhline(0,ls=":")
axvline(0,ls=":")
grid()
xlabel(r"$\theta_{obs}$")
ylabel("95% limit")
Out[20]:
Text(0,0.5,'95% limit')

Situace změní, pokud v Bayesovském přístupu zvolíme jinou a-priorní pravděpodobnost $\pi(\theta)$, např. podle Jeffrey $1/x$...

In [21]:
at0=0.1
jeff_gauss=lambda x,m:exp(-(x-m)**2)/x
jnorm=r_[[integrate.quad(jeff_gauss,at0,20,i)[0] for i in cpos]]
plot(cpos,jnorm)
xlabel("skut. hodnota")
ylabel("integr. 0-infty (Jeffrey prior)")
Out[21]:
Text(0,0.5,'integr. 0-infty (Jeffrey prior)')
In [23]:
w1=lambda p,x:p[0]*exp(-(x-p[1])**2/2./p[2]**2)
w3=lambda p,x:p[0]*exp((-(x-p[1])**2+p[3]*(x-p[1])**3)/2./p[2]**2)
normmod=optimize.leastsq(lambda p:w1(p,cpos)-jnorm,[jnorm.sum(),0,1.])[0]
normmod3=optimize.leastsq(lambda p:w3(p,cpos)-jnorm,[jnorm.sum(),0,1.,0])[0]
plot(cpos,w1(normmod,cpos)-jnorm)
plot(cpos,w3(normmod3,cpos)-jnorm)
grid()
title("aproximace")
ylabel("reziduum")
legend(["w1: gaussovka","w3: i clen 3.radu"])
Out[23]:
<matplotlib.legend.Legend at 0x1ff8e952240>

vzhledem k asymetrii neni gausovka vhodna funkce

In [33]:
jvars=[[integrate.quad(jeff_gauss,at0,k,i)[0] for i in cpos] for k in r_[0.4:2.5:.4]]
#ruzne horni meze 0.4-2.5
[plot(cpos,b) for b in jvars]
plot(cpos,jnorm,':')
xlabel("skut. hodnota")
ylabel("integ 0.1 -- horni mez")
legend(r_[0.4:2.5:.4],loc=2)
Out[33]:
<matplotlib.legend.Legend at 0x1ff8ed35dd8>

závislost parametru gaussovky na horní mezi se pokoušíme modelovat kvadratickými polynomy

In [27]:
# nadale se snazime fitovat gausovky
varep=r_[[optimize.leastsq(lambda p:w3(p,cpos)-array(b),[sum(b),0,1.,0])[0] for b in jvars]]
fullmod=[polyfit(r_[0.4:2.5:.4],v,2) for v in r_[varep].transpose()]
fullmod
Out[27]:
[array([-0.45439472,  1.75532499,  0.80279391]),
 array([-0.08246034,  0.38448172,  0.07608313]),
 array([0.00816365, 0.06968679, 0.67869994]),
 array([ 0.03809857, -0.01861811, -0.00059378])]
In [28]:
# jak dobry je fit 
[(abs(polyval(fullmod[i],r_[0.4:2.5:.4])-varep[:,i])).max() for i in range(4)]
Out[28]:
[0.08746860198906292,
 0.0020304554895774762,
 0.006849261007292995,
 0.006491470118218573]
In [36]:
i=3
plot(r_[0.4:2.5:.4],polyval(fullmod[i],r_[0.4:2.5:.4]))
plot(r_[0.4:2.5:.4],varep[:,i],'+')
ylabel("asym. clen")
title("aproximace parametru fitu asym. gaussovkou")
Out[36]:
Text(0.5,1,'aproximace parametru fitu asym. gaussovkou')
In [42]:
# horni limit v interv. [0.4:2.5]
# param. p(stred rozdeleni) v interv. [-3:2]
confid=lambda upl,cent:w3([polyval(f,upl) for f in fullmod],cent)/w3(normmod3,cent)
#funkce modelující poměr integrálů (věrohodností)
In [21]:
#poměr vůči integrálu přes kladná čísla pro ruzne hodnoty stredu rozdeleni (a uplim 2.0)
jvars[4]/jnorm
Out[21]:
array([ 0.99999999,  0.99999998,  0.99999997,  0.99999996,  0.99999994,
        0.99999992,  0.99999988,  0.99999983,  0.99999975,  0.99999963,
        0.99999946,  0.99999921,  0.99999885,  0.99999832,  0.99999755,
        0.99999642,  0.99999479,  0.99999241,  0.99998895,  0.99998391,
        0.99997661,  0.999966  ,  0.99995063,  0.99992836,  0.99989615,
        0.99984962,  0.99978253,  0.99968594,  0.99954717,  0.99934829,
        0.99906404,  0.99865901,  0.99808399,  0.997271  ,  0.99612708,
        0.99452641,  0.99230091,  0.98922946,  0.98502638,  0.9793306 ,
        0.9716977 ,  0.96159794,  0.94842445,  0.93151587,  0.91019708,
        0.88383912,  0.85193529,  0.81418517,  0.77057294,  0.72142424])
In [52]:
#bayes solution
uplims=[.5]
lolims=[4.]
for i in r_[-2.5:4:0.5]:
    lolims.append(optimize.fsolve(lambda p:confid(p,i)-0.1,lolims[-1])[0])
    uplims.append(optimize.fsolve(lambda p:confid(p,i)-0.9,uplims[-1])[0])# 90 %
del uplims[0]
del lolims[0]
In [54]:
plot(r_[-2.5:4:0.5],uplims)
plot(r_[-2.5:4:0.5],lolims)
grid()
In [70]:
xbas=r_[0.1:2.5:.2]
plot(xbas,[confid(q,3) for q in xbas],hold=0)
plot(xbas,[confid(q,2) for q in xbas],hold=1)
plot(xbas,[confid(q,1) for q in xbas],hold=1)
plot(xbas,[confid(q,0.) for q in xbas],hold=1)
grid()
# flat (nonnegative) prior
plot(xbas,[1-stats.norm(3,1).sf(q)/stats.norm(3,1).sf(0) for q in xbas],'b:')
plot(xbas,[1-stats.norm(2,1).sf(q)/stats.norm(2,1).sf(0) for q in xbas],'g:')
plot(xbas,[1-stats.norm(1,1).sf(q)/stats.norm(1,1).sf(0) for q in xbas],'r:')
plot(xbas,[1-stats.norm(0,1).sf(q)/stats.norm(0,1).sf(0) for q in xbas],'c:')

title("jeffrey's prior vs. flat")
legend(["cent=3","cent=2","cent=1","cent=0"],loc=4)
xlabel("upper limit")
ylabel("frac to total")
ylim(-.2,1.1)
Out[70]:
(-0.2, 1.1)
In [71]:
#prime nalezeni horni meze 
zept=lambda q,lp:1-stats.norm(2,1).sf(q)/stats.norm(2,1).sf(0)-lp
optimize.fsolve(zept,3,0.9)
Out[71]:
array([ 3.29462399])

Horni / dolni limita s rovnom. apriorni pravdepodobnosti

In [73]:
xbas=r_[-2.5:4:0.5]
lolims2=r_[[stats.norm(i,1).isf(stats.norm(i,1).sf(0)*0.9) for i in xbas]]
uplims2=r_[[stats.norm(i,1).isf(stats.norm(i,1).sf(0)*0.1) for i in xbas]]
plot(xbas,uplims2,hold=0)
plot(xbas,lolims2,hold=1)
Out[73]:
[<matplotlib.lines.Line2D at 0xe2328d0>]
In [3]:
import sys;sys.path.append("/Users/Admin/Code/")
import nbexport as ne
ne.odir="/Users/Admin/Documents/Vyuka/MMZM/export/"
ne.export("Index.ipynb")

saved /Users/Admin/Documents/Vyuka/MMZM/export/Index.html+ 0 image(s)
ok
Out[3]:
51