Horní limity (pokračování)

pro normální rozdělení, a-priorní bayesovská funkce podle Cauchyho (nemá problém v 0, integrál konečný) pro kladná $\nu$:

$$\int_0^{\nu^{up}}{P(\nu|x_{obs})} d \nu = \frac{\int_0^{\nu^{up}} L(x_{obs}|\nu) \pi(\nu) \ d \nu}{\int_0^{\infty} L(x_{obs}|\nu) \pi(\nu) \ d \nu}= \frac{\int_0^{\nu^{up}} \exp(-(x_{obs}-\nu)^2) \frac{1}{1+\nu^2} \ d \nu}{\int_0^{\infty} \exp(-(x_{obs}-\nu)^2) \frac{1}{1+\nu^2} \ d \nu}$$

In [2]:
from numpy import *
from matplotlib.pyplot import *
%matplotlib inline
from scipy import stats,optimize,integrate
cauc_gauss=lambda x,m:exp(-(x-m)**2)/(x**2+1)
at0=0
cpos=r_[-2:4:0.1]
# vypocet integralu 
jnorm=r_[[integrate.quad(cauc_gauss,at0,20,i)[0] for i in cpos]]
In [3]:
plot(cpos,jnorm)
xlabel('merene x')
ylabel('1-beta')
Out[3]:
<matplotlib.text.Text at 0x8e701d0>

aproximujeme funkci pro další použití Gram-Charlierovou modifikací Gaussovy funkce

$$\exp\frac{(x-\mu)^2}{2\sigma^2}\ \left( 1+\sum_i {\alpha_i H_i(x)} \right)$$

kde H(i) jsou Hermiteovy polynomy (ortogonální s příslušnou vahou)

In [4]:
a0=[2,3,1,1,1,1.01]
from scipy import special
#gram-charlier modification
gram=lambda x,a:(sum([special.hermite(i+1)*a[i] for i in range(len(a)-4)])+1)(x)*a[-1]*exp(-(x-a[-3])**2/2/a[-2]**2)
In [5]:
a=optimize.leastsq(lambda a:gram(cpos,a)-jnorm,a0)[0]
In [6]:
plot(cpos,gram(cpos,a))
grid()
In [7]:
# rozdil obou funkci
a0=[ 0.01311388, -0.06893304,  0.01198881, 0.,  1.15141975,  1.094612  ,  1.01428982]
plot(cpos,gram(cpos,a0)-jnorm)
grid()
a
Out[7]:
array([-0.33577522,  0.03837713,  1.        ,  1.52678409,  0.97090018,
        2.43587294])
In [10]:
jvars=[[integrate.quad(cauc_gauss,at0,k,i)[0] for i in cpos] for k in r_[0.4:2.5:.4]]
#integraly pro ruzne hodnoty horni meze (0.4 az 2.5)
[plot(cpos,b) for b in jvars]
legend(r_[0.4:2.5:.4])
Out[10]:
<matplotlib.legend.Legend at 0x9e1d080>
In [23]:
rep=[[.02,-.02,.5,1.,1.,1.]]
sum(gram(cpos,rep[-1]))
Out[23]:
23.055798673680464
In [24]:
#varep2=r_[[optimize.leastsq(lambda a:gram(cpos,a)-array(b),a0)[0] for b in jvars]]
rep=[[.02,-.02,0,.5,1.,1.,1.]]
#rep=[[.3,0,1.]]
dis=[]
for i in range(len(jvars)):
#    rfit=optimize.leastsq(lambda a:a[0]*exp(-(cpos-a[1])**2/2/a[2]**2)-array(jvars[i]),rep[-1])
    rfit=optimize.leastsq(lambda a:gram(cpos,a)-array(jvars[i]),rep[-1])
    rep.append(rfit[0])
    a=rfit[0]
    model=gram(cpos,a)#a[0]*exp(-(cpos-a[1])**2/2/a[2]**2)
    dis.append([abs(model-jvars[i]).max(),((model-jvars[i])**2).sum()])
In [28]:
rep=array(rep)
rep[:,-1]
#(sum([special.hermite(i+1)*a[i] for i in range(len(a)-4)])+1)(3)*exp(-(3-a[-3])**2/2/a[-2]**2
a,rep[1]
Out[28]:
(array([ 38.87939625,   0.77865293,  -0.83484855,   0.5       ,
        -83.30275575,   1.25004779, -82.45115599]),
 array([ 38.87939625,   0.77865293,  -0.83484855,   0.5       ,
        -83.30275575,   1.25004779, -82.45115599]))

Poissonovská data s pozadím

$L(\nu_s)=(\nu_s+\nu_b)^{n_{obs}} \exp(-\nu_s-\nu_b)$, kde $\nu_b$ je pozadí

Bayesovsky (rovnoměrná a-priorní fce pro x>0): $$1-\beta = \frac{\int_0^{\nu_s^{up}} L(n_{obs}|\nu_s)\ d \nu_s}{\int_0^{\infty} L(n_{obs}|\nu_s)\ d \nu_s}= \frac{\int_0^{\nu_s^{up}} (\nu_s+\nu_b)^{n_{obs}} \exp(-\nu_s-\nu_b)\ d \nu_s}{\int_0^{\infty} (\nu_s+\nu_b)^{n_{obs}} \exp(-\nu_s-\nu_b)\ d \nu_s} $$

Dolní integrál lze spočítat per-partes ($n_{obs}$ krát) jako sumu

$$\exp(-\nu_b) \sum_{n=0}^{n_{obs}}\frac{(\nu_b)^n}{n!}$$ místo horního vzít integrál od $\nu_s^{up}$ do $\infty$ a dostaneme

$$\beta=\frac{\exp(-\nu_s^{up}-\nu_b) \sum_{n=0}^{n_{obs}} (\nu_s^{up}+\nu_b)^n / n!}{\exp(-\nu_b) \sum_{n=0}^{n_{obs}} (\nu_b)^n / n!}$$

Existuje i matem. vztah pro uvednou řadu využívající kumul. distribuční funkci $\chi^2$ rozdělení

$$ \exp(-\nu) \sum_{n=0}^{n_{obs}} \frac{\nu^n}{n!} = \int_{2\nu}^{\infty} f_{\chi^2}(z;d.o.f.=2(n_{obs}+1))\ dz = 1-F_{\chi^2}(2\nu;d.o.f.=2(n_{obs}+1))$$

In [16]:
from scipy.misc import factorial

betsum=lambda nup,nobs:exp(-nup)*(pow(nup,r_[:nobs+1])/factorial(r_[:nobs+1])).sum()
nub=2
aobs=3
norm=betsum(nub,aobs)
betas=lambda nup:betsum(nup+nub,aobs)/norm
In [17]:
nux=r_[2:5:.1]
plot(nux,[betas(x) for x in nux])
aobs=4
norm=betsum(nub,aobs)
betas=lambda nup:betsum(nup+nub,aobs)/norm
plot(nux,[betas(x) for x in nux])
title('pozadi=2')
xlabel('horni limita')
ylabel('beta')
legend(['nobs=3','nobs=4'])
Out[17]:
<matplotlib.legend.Legend at 0xa57f940>
In [18]:
from scipy import optimize
def get_uplim(beta=0.1,nub=2,nobs=3,ninit=6):
    norm=betsum(nub,nobs)
    betas=lambda nup:betsum(nup+nub,nobs)/norm
    return optimize.fsolve(lambda x:betas(x)-beta,ninit)[0]
rep=[get_uplim(nobs=4,nub=i) for i in range(1,9)]
rep
Out[18]:
[6.99997179468492,
 6.0874363960126416,
 5.3446841011119481,
 4.7792510521283305,
 4.3570003069452019,
 4.0399668739655832,
 3.7981024269423505,
 3.6099668723410896]
In [10]:
idx=polyfit(r_[1:9],log(r_[rep]-3),1)
plot(r_[1:9],rep)
plot(r_[1:9],exp(polyval(idx,r_[1:9]))+3)
Out[10]:
[<matplotlib.lines.Line2D at 0x2f78710>]