누적 확률 (Python)이 주어지면 랜덤 변수의 값을 얻습니다.

Question

다음은 간단한 배경 정보입니다. Monte-Carlo 접근법을 사용하여 두 로그 정규 랜덤 변수의 선형 조합에 대해 결합 된 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. 컨볼 루션 과 유사하게 푸리에 변환을 앞뒤로 수행

4. 물론 Monte Carlo 통합은 항상 옵션입니다.

최신 정보

여기에 가장 간단한 코드가 있습니다.

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

그래프

Answer 2

2 로그 정규 분포의 합에서 표본을 얻으려면 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]))

x1과 x2의 합은 그 자체가 분포입니다.

sum = x1+x2

평균 sum.getMean()[0](= 1.5379) 또는 표준 편차 sum.getStandardDeviation()[0](= 0.42689241033309544)에 액세스 할 수 있습니다.

물론 N = 5에 대해 N 크기의 샘플을 얻을 수 있습니다. sum.getSample(5)

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