संचयी संभाव्यता (पायथन) को दिए गए यादृच्छिक चर के मूल्य को प्राप्त करें

Aug 16 2020

यहाँ एक त्वरित पृष्ठभूमि जानकारी है। मैं मोंटे-कार्लो दृष्टिकोण का उपयोग करते हुए दो lognormal यादृच्छिक चर के रैखिक संयोजन के लिए एक संयुक्त सीडीएफ प्राप्त करने की कोशिश कर रहा हूं और फिर, नमूना करने के लिए इसे उल्टा कर देता हूं। यहाँ वही करने के लिए पायथन कोड दिया गया है:

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

प्रशन:

क्या इस ऑपरेशन को करने के लिए घबराहट में एक सीधा कार्य है?

कोड की अंतिम पंक्ति में, मैं माध्य मान ले रहा हूं, क्या कोई तरीका है जिससे मैं प्रक्षेप, आदि द्वारा अधिक सटीक मान प्राप्त कर सकता हूं? यदि हां, तो मैं इसे पायथन में कैसे लागू कर सकता हूं

जवाब

3 SeverinPappadeux Aug 16 2020 at 23:36

जब आप दो आरवी एक्स + वाई, पीडीएफ एक्स (एक्स), पीडीएफ वाई (वाई) जानते हैं और पीडीएफ एक्स + वाई (जेड) जानना चाहते हैं , तो अच्छी तरह से ज्ञात मामला है । आप यहां समान दृष्टिकोण का उपयोग कर सकते हैं, पीडीएफ की गणना कर सकते हैं और सीडीएफ = डी पीडीएफ (जेड) / डीजेड बना सकते हैं

PDF aX + bY (z) = S डाई पीडीएफ Y (y) पीडीएफ X ((z-by) / a) / a |

जहां Sएकीकरण को दर्शाता है।

आप इसे सीधे CDF के लिए लिख सकते हैं

CDF aX + bY (z) = S डाई पीडीएफ Y (y) CDF X ((z-by) / a)

आप इस अभिन्न की गणना कर सकते हैं:

  1. विश्लेषणात्मक

  2. संख्यात्मक रूप से, SciPy का उपयोग करते हुए

  3. कॉन्ट्रोवर्सी के समान ही फूरियर आगे और पीछे ट्रांसफॉर्म होता है

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

ग्राफ़

JeanA. Oct 22 2020 at 00:55

यदि आप 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]))

एक्स 1 और एक्स 2 का योग अपने आप में एक वितरण है

sum = x1+x2

आप इसके माध्य sum.getMean()[0](= 1.5379) या इसके मानक विचलन sum.getStandardDeviation()[0](= 0.42689241033309544) पर पहुँच सकते हैं

और निश्चित रूप से, आप N = 5 के लिए किसी भी आकार का एक नमूना प्राप्त कर सकते हैं: sum.getSample(5)

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