Dapatkan nilai variabel acak yang diberi probabilitas kumulatif (Python)

Aug 16 2020

Berikut adalah informasi latar belakang singkat. Saya mencoba untuk mendapatkan CDF gabungan untuk kombinasi linier dari dua variabel acak lognormal menggunakan pendekatan Monte-Carlo dan kemudian, balikkan untuk melakukan pengambilan sampel. Berikut kode Python untuk melakukan hal yang sama:

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])

Pertanyaan:

Apakah ada fungsi langsung di scipy untuk melakukan operasi ini?

Di baris terakhir kode, saya mengambil nilai rata-rata, apakah ada cara saya bisa mendapatkan nilai yang lebih akurat dengan interpolasi, dll? Jika demikian, Bagaimana cara menerapkannya dengan Python

Jawaban

3 SeverinPappadeux Aug 16 2020 at 23:36

Nah, ada kasus terkenal ketika Anda menjumlahkan dua RV X + Y, tahu PDF X (x), PDF Y (y) dan ingin tahu PDF X + Y (z). Anda dapat menggunakan pendekatan serupa di sini, menghitung PDF dan membuat CDF = d PDF (z) / dz

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

dimana Smenunjukkan integrasi.

Anda dapat menuliskannya langsung untuk CDF

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

Anda dapat menghitung integral ini:

  1. Secara analitis

  2. Secara numerik, menggunakan SciPy

  3. Lakukan transformasi Fourier maju dan mundur, mirip dengan Konvolusi

  4. Tentu saja, integrasi Monte Carlo selalu menjadi pilihan

MEMPERBARUI

Berikut adalah kode paling sederhana untuk membantu Anda memulai

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

Jika Anda ingin mendapatkan sampel dari jumlah 2 distribusi lognormal, Anda tidak memerlukan skema Monte-Carlo.

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]))

jumlah dari x1 dan x2 itu sendiri adalah sebuah distribusi

sum = x1+x2

Anda dapat mengakses mean sum.getMean()[0](= 1,5379) atau deviasi standarnya sum.getStandardDeviation()[0](= 0,42689241033309544)

dan tentu saja, Anda bisa mendapatkan sampel dengan ukuran berapa pun N Untuk N = 5: sum.getSample(5)

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