Kombinace náhodných proměnných

Simulace součinu dvou normálně rozdělených veličin
numerický výpočet konvoluce (Gaussova kvadratura)

  • interval $(-\infty,\infty)$ nahrazen $(-100,100)$
In [1]:
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
In [2]:
x=r_[0.1:5:0.1]
fun=multiple()
plot(x,[fun(p) for p in x])
Out[2]:
[<matplotlib.lines.Line2D at 0x3607250>]

analytický výsledek

modifikovaná Besselova funkce 2. druhu ($K_0$)

In [9]:
from scipy import special
x=r_[0.05:5:.05]
plot(x,special.k0(abs(x)))
#special.k0?
In [4]:
v1=normal(size=500)
v12=normal(size=500)
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í součinu necentrálních Gaussovek

posun vpravo o 1,2,3

In [5]:
w3a=hist((v1+1)*(v12+1),p,alpha=0.4)
w3b=hist((v1+2)*(v12+2),p,alpha=0.4)
w3c=hist((v1+3)*(v12+3),p,alpha=0.4)

Výpočet šikmosti a špičatosti pro součin 2 norm. rozd. veličin s postupně rostoucí střed. hodnotou

In [6]:
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(1,11)]
array(comps).transpose()
Out[6]:
array([[  1.05848776e+00,   4.02813798e+00,   8.99778820e+00,
          1.59674384e+01,   2.49370886e+01,   3.59067389e+01,
          4.88763891e+01,   6.38460393e+01,   8.08156895e+01,
          9.97853397e+01],
       [  2.86115060e+00,   9.21792290e+00,   2.00622797e+01,
          3.53942209e+01,   5.52137467e+01,   7.95208569e+01,
          1.08315552e+02,   1.41597831e+02,   1.79367694e+02,
          2.21625142e+02],
       [  1.11334344e+00,   7.45462715e-01,   4.85368096e-01,
          3.34606934e-01,   2.39364178e-01,   1.74391506e-01,
          1.27432433e-01,   9.19806382e-02,   6.42995981e-02,
          4.21019188e-02],
       [  2.78118917e+00,   7.81247918e-01,   2.27421146e-01,
          2.74663286e-02,  -6.24346271e-02,  -1.09044457e-01,
         -1.35660647e-01,  -1.51940176e-01,  -1.62414239e-01,
         -1.69415512e-01]])
In [7]:
acomp=array(comps).transpose()
a=range(1,11)
plot(a,sqrt(acomp[1])/a)
plot(a,acomp[2])
plot(a,acomp[3])
legend(["disp./a","skew","kurt"])
xlabel("necentralnost")
grid()

"Non-central" Chi square

nahodna promenna: (gaus + A)*(gaus + A)=A^2 + 2A gaus + chi2(1 d.f.)

$\exp(-(x-A)^2/2/s^2)=\exp(-A^2/2/s^2)*\exp(xA/s^2)*\exp(-x^2/2/s^2)$

In [10]:
v2=normal(loc=10,size=1000)
v4=normal(loc=14,size=1000)

w1=hist(v2,10)
w2=hist(v2**2,20)
w4=hist(v4**2,20)