รับค่าของตัวแปรสุ่มตามความน่าจะเป็นสะสม (Python)
นี่คือข้อมูลพื้นฐานโดยย่อ ฉันกำลังพยายามหา 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 ได้อย่างไร
คำตอบ
มีกรณีที่รู้จักกันดีเมื่อคุณรวม RVs 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)
คุณสามารถคำนวณอินทิกรัลนี้:
ในเชิงวิเคราะห์
ในเชิงตัวเลขโดยใช้SciPy
ทำการแปลงฟูเรียร์ไปข้างหน้าและข้างหลังคล้ายกับConvolution
แน่นอนว่าการรวม 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()
กราฟ

หากคุณต้องการรับตัวอย่างจากผลรวมของการแจกแจงแบบ lognormal 2 ตัวคุณไม่จำเป็นต้องมีโครงร่างมอนติคาร์โล
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 สำหรับ N = 5: sum.getSample(5)
print(sum.getSample(5))
0 : [ 1.29895 ]
1 : [ 1.32224 ]
2 : [ 1.259 ]
3 : [ 1.16083 ]
4 : [ 1.30129 ]