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}$$
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]]
plot(cpos,jnorm)
xlabel('merene x')
ylabel('1-beta')
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)
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)
a=optimize.leastsq(lambda a:gram(cpos,a)-jnorm,a0)[0]
plot(cpos,gram(cpos,a))
grid()
# 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
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])
rep=[[.02,-.02,.5,1.,1.,1.]]
sum(gram(cpos,rep[-1]))
#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()])
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]
$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))$$
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
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'])
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
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)