Kümülatif olasılık (Python) verildiğinde rastgele değişkenin değerini elde edin

Aug 16 2020

İşte hızlı bir arka plan bilgisi. Monte-Carlo yaklaşımını kullanarak iki lognormal rastgele değişkenin doğrusal kombinasyonu için birleşik bir CDF elde etmeye ve ardından örnekleme yapmak için onu tersine çevirmeye çalışıyorum. İşte aynısını yapmak için Python kodu:

import numpy as np
from scipy import special


# parameters of distribution 1
mu1 = 0.3108
s1=0.3588

# parameters of distribution 2
mu2=1.2271
s2=0.2313

a = 2
b=3

N_sampling = 10000

kk=0

Y=np.zeros(N_sampling)
X1=np.zeros(N_sampling)
X2=np.zeros(N_sampling)

while(kk<N_sampling):
    F = np.random.rand(2)
    X1[kk]=np.exp(mu1+(2**0.5)*s1*special.erfinv(2*F[0]-1))  # sampling X1 (distribution1) by inverting the CDF
    X2[kk]=np.exp(mu2+(2**0.5)*s2*special.erfinv(2*F[1]-1))  # sampling X2 (distribution2) by inverting the CDF  
    
    Y[kk]=a*X1[kk]+b*X2[kk] # obtain the random variable as a linear combination of X1 and X2
    kk=kk+1
    

# Obtain the CDF of Y

freq, bin_borders = np.histogram(Y, bins=50)    
norm_freq = freq/np.sum(freq)
cdf_Y = np.cumsum(norm_freq)


# obtain the value of Y given the value of cdf_Y
cdf_Y_input=0.5
idx=np.searchsorted(cdf_Y,cdf_Y_input)
Y_out = 0.5*(bin_borders[idx-1]+bin_borders[idx])

Sorular:

Scipy'de bu işlemi gerçekleştirmek için doğrudan bir işlev var mı?

Kodun son satırında ortalama değeri alıyorum, enterpolasyon ile daha doğru değerler elde etmenin bir yolu var mı? Öyleyse, Python'da nasıl uygularım

Yanıtlar

3 SeverinPappadeux Aug 16 2020 at 23:36

İki RV'yi X + Y topladığınızda, PDF X (x), PDF Y (y) bildiğiniz ve PDF X + Y (z) bilmek istediğinizde iyi bilinen bir durum vardır . Burada benzer bir yaklaşım kullanabilir, PDF'yi hesaplayabilir ve CDF = d PDF (z) / dz yapabilirsiniz.

PDF aX + bY (z) = S dy PDF Y (y) PDF X ((z-by) / a) / | a |

burada Sentegre belirtmektedir.

Doğrudan CDF için yazabilirsiniz

CDF aX + bY (z) = S dy PDF Y (y) CDF X ((z-by) / a)

Bu integrali hesaplayabilirsiniz:

  1. Analitik olarak

  2. Sayısal olarak, SciPy kullanarak

  3. Fourier, Convolution'a benzer şekilde ileri ve geri dönüşüm mü?

  4. Elbette, Monte Carlo entegrasyonu her zaman bir seçenektir

GÜNCELLEME

İşte sizi harekete geçirecek en basit kod

import numpy as np
from math import erf

SQRT2 = np.sqrt(2.0)
SQRT2PI = np.sqrt(2.0*np.pi)
    
def PDF(x):
    if x <= 0.0:
        return 0.0

    q = np.log(x)
    return np.exp( - 0.5*q*q ) / (x * SQRT2PI)

def CDF(x):
    if x <= 0.0:
        return 0.0

    return 0.5 + 0.5*erf(np.log(x)/SQRT2)    

import scipy.integrate as integrate
import matplotlib.pyplot as plt

a = 0.4
b = 0.6

N = 101

z = np.linspace(0.0, 5.0, N)
c = np.zeros(N) # CDF of the sum
p = np.zeros(N) # PDF of the sum
t = np.zeros(N) # CDF as integral of PDF

for k in range(1, N):
    zz = z[k]
    ylo = 0.0
    yhi = zz/b

    result = integrate.quad(lambda y: PDF(y) * CDF((zz - b*y)/a), ylo, yhi)
    print(result)
    c[k] = result[0]

    result = integrate.quad(lambda y: PDF(y) * PDF((zz - b*y)/a)/a, ylo, yhi)
    print(result)
    p[k] = result[0]

    t[k] = integrate.trapz(p, z) # trapezoidal integration over PDF


plt.plot(z, c, 'b^') # CDF
plt.plot(z, p, 'r.') # PDF
plt.plot(z, t, 'g-') # CDF as integral over PDF
plt.show()

Grafik

JeanA. Oct 22 2020 at 00:55

2 lognormal dağılımın toplamından bir örnek almak istiyorsanız, bir Monte-Carlo şemasına ihtiyacınız yoktur.

import openturns as ot 
x1 = ot.LogNormal()
x1.setParameter(ot.LogNormalMuSigma()([0.3108, 0.3588, 0.0]))
# in order to convert mu, sigma into mulog and sigmalog

x2 = ot.LogNormal()
x2.setParameter(ot.LogNormalMuSigma()([1.2271, 0.2313, 0.0]))

x1 ve x2'nin toplamının kendisi bir dağılımdır

sum = x1+x2

ortalamasına sum.getMean()[0](= 1.5379) veya standart sapmasına sum.getStandardDeviation()[0](= 0.42689241033309544) erişebilirsiniz

ve tabii ki, N = 5 için herhangi bir N boyutundan bir örnek alabilirsiniz: sum.getSample(5)

print(sum.getSample(5))
0 : [ 1.29895 ]
1 : [ 1.32224 ]
2 : [ 1.259   ]
3 : [ 1.16083 ]
4 : [ 1.30129 ]