python - 获取给定累积概率的随机变量的值 (Python)

Question

这是一个快速的背景信息。我试图使用蒙特卡罗方法获得两个对数正态随机变量的线性组合的组合 CDF，然后将其反转以进行采样。下面是执行相同操作的 Python 代码:

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

问题:

scipy中有直接函数来执行此操作吗？

在代码的最后一行，我取平均值，有没有办法通过插值等获得更准确的值？如果是这样，我如何在Python中实现它

Answer 1

嗯，有一个众所周知的情况，当您将两个 RV X+Y 相加，知道 PDF_X(x)、PDF_Y(y) 并且想知道PDF_X+Y(z)。您可以在这里使用类似的方法，计算 PDF 并使 CDF=d PDF(z)/dz

PDF_aX+bY(z) = S dy PDF_Y(y) PDF_X((z-by)/a)/|a|

其中S表示集成。

可以直接写成CDF

CDF_aX+bY(z) = S dy PDF_Y(y) CDF_X((z-by)/a)

你可以计算这个积分:

1. 分析

2. 数字上，使用 SciPy

3. 向前和向后进行傅里叶变换，类似于 Convolution

4. 当然，蒙特卡罗积分始终是一种选择

更新

这是最简单的代码，可以帮助您开始工作

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

图表