In [3]:
x=r_[0.1:5:0.1]
fun=multiple()
prof=[fun(r) for r in x]
plot(x,prof)
Out[3]:
Simulace podílu dvou normálně rozdělených veličin
numerický výpočet konvoluce (Gaussova kvadratura)
%pylab inline
from scipy import integrate
def multiple(wid=3,infty=100):
sgau=lambda x:exp(-x**2/wid**2)
n1=integrate.quad(sgau,-infty,infty)
ssqr=lambda x,z: sgau(z*x)*sgau(x)*abs(x)
smul=lambda z:integrate.quad(ssqr,-infty,infty+10,(z,))[0]/n1[0]**2
return smul
x=r_[0.1:5:0.1]
fun=multiple()
prof=[fun(r) for r in x]
plot(x,prof)
v případě centrovaných ($E[x]=0$) norm. rozdělení je podílem Cauchyho rozdělení
from scipy import special
#x=r_[0.1:5:.1]
plot(x,prof-1/(1+abs(x)**2)/pi)
ylabel("rozdil numericky-analyticky")
#special.k0?
x=r_[-5:5:0.05]
[plot(x,g/(g**2+(x-1)**2)/pi) for g in [0.7,1,2]]
grid()
v1=normal(size=1000)
v12=normal(size=1000)
p=r_[-10:10:0.4]
w1=hist(v1,p,alpha=0.3)
w2=hist(v1**2,p,alpha=0.3)
w3=hist(v1*v12,p,alpha=0.3)
porovnání podílu necentrálních Gaussovek
posun vpravo o 10,15,20
p=r_[0.5:1.5:0.02]
w3a=hist((v1+10)/(v12+10),p,alpha=0.4)
w3b=hist((v1+15)/(v12+15),p,alpha=0.4)
w3c=hist((v1+20)/(v12+20),p,alpha=0.4)
aproximace z Gaussovy křivky
mx,my=8,8
sx,sy=1,1
vG=lambda z:(my-mx*z)/sqrt(sx**2*z**2+sy**2)
bG=lambda z:(mx*sy/sx+my*sx/sy*z)/sqrt(sx**2*z**2+sy**2)
pZ=lambda z,b:sx*sy/(sx**2*z**2+sy**2)*(1+b*exp(b**2/2)*(special.erf(b/sqrt(2))-0.5)*sqrt(pi/2))*exp(-(mx**2/sx**2+my**2/sy**2)/2.)/pi**2/2.
x=r_[0.1:2:.02]
plot(x,exp(-vG(x)**2/2)/sqrt(2*pi))
mx,my=15,8
plot(x,exp(-vG(x)**2/2)/sqrt(2*pi))
plot(x,pZ(x,bG(x)))
def moms(xd):
xm=xd.mean()
x2=((xd-xm)**2).mean()
x3=((xd-xm)**3).mean()/pow(x2,1.5)
x4=((xd-xm)**4).mean()/x2**2-3
return xm,x2,x3,x4
comps=[moms((v1+i)/(v12+i)) for i in range(5,15)]
print("2.moment",array(comps)[:,0])
print("3.moment",array(comps)[:,1])
print("4.moment",array(comps)[:,1])