Nhận giá trị của biến ngẫu nhiên theo xác suất tích lũy (Python)

Question

Đây là một thông tin cơ bản nhanh chóng. Tôi đang cố gắng lấy CDF kết hợp cho sự kết hợp tuyến tính của hai biến ngẫu nhiên bất thường bằng cách sử dụng phương pháp tiếp cận Monte-Carlo và sau đó, đảo ngược nó để lấy mẫu. Đây là mã Python để làm điều tương tự:

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

Câu hỏi:

Có một chức năng trực tiếp trong scipy để thực hiện thao tác này không?

Trong dòng cuối cùng của mã, tôi đang lấy giá trị trung bình, có cách nào tôi có thể nhận được giá trị chính xác hơn bằng cách nội suy, v.v. không? Nếu vậy, làm cách nào để triển khai nó bằng Python

Answer 1

Vâng, có một trường hợp nổi tiếng khi bạn tính tổng hai RV X + Y, biết PDF _X (x), PDF _Y (y) và muốn biết PDF _{X + Y} (z). Bạn có thể sử dụng phương pháp tương tự ở đây, tính toán PDF và tạo CDF = d PDF (z) / dz

PDF _{aX + bY} (z) = S dy PDF _Y (y) PDF _X ((z-by) / a) / | a |

nơi Sbiểu thị sự tích hợp.

Bạn có thể viết nó trực tiếp cho CDF

CDF _{aX + bY} (z) = S dy PDF _Y (y) CDF _X ((z-by) / a)

Bạn có thể tính tích phân này:

1. Phân tích

2. Về số lượng, sử dụng SciPy

3. Do Fourier biến đổi tiến và lùi, tương tự như Convolution

4. Tất nhiên, tích hợp Monte Carlo luôn là một lựa chọn

CẬP NHẬT

Đây là mã đơn giản nhất để giúp bạn tiếp tục

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

Đồ thị

Answer 2

Nếu bạn muốn lấy một mẫu từ tổng của 2 phân phối chuẩn, bạn không cần lược đồ 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]))

tổng của x1 và x2 chính nó là một phân phối

sum = x1+x2

bạn có thể truy cập giá trị trung bình sum.getMean()[0](= 1.5379) hoặc độ lệch chuẩn của nó sum.getStandardDeviation()[0](= 0.42689241033309544)

và tất nhiên, bạn có thể lấy một mẫu có kích thước N bất kỳ Với N = 5: sum.getSample(5)

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