Obtain the value of the random variable given the cumulative probability (Python)

user_5129061 0
Viewed 214

Here is a quick background information. I am trying to obtain a combined CDF for the linear combination of two lognormal random variables using Monte-Carlo approach and then, invert it to do sampling. Here is the Python code for doing the same:

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

Questions:

Is there a direct function in scipy to perform this operation ?

In the last line of the code, I am taking the mean value, is there a way I could obtain more accurate values by interpolation, etc ? If so, How do I implement it in Python

probability-distribution probability distribution statistics python
2 Answers

Well, there is well-known case when you sum two RVs X+Y, know PDFX(x), PDFY(y) and want to know PDFX+Y(z). You could use similar approach here, compute PDF and make CDF=d PDF(z)/dz

PDFaX+bY(z) = S dy PDFY(y) PDFX((z-by)/a) / |a|

where S denotes integration.

You could write it directly for CDF

CDFaX+bY(z) = S dy PDFY(y) CDFX((z-by)/a)

You could compute this integral:

  1. Analytically

  2. Numerically, using SciPy

  3. Do Fourier transform forward and backward, similar to Convolution

  4. Of course, Monte Carlo integration is always an option

UPDATE

Here is the simplest code to get you going

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

Graph

enter image description here

user_40446960

If you want to get a sample from the sum of 2 lognormal distribution, you do not need a Monte-Carlo scheme.

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

the sum of x1 and x2 is itself a distribution

sum = x1+x2

you can access its mean sum.getMean()[0] (= 1.5379) or its standard deviation sum.getStandardDeviation()[0](= 0.42689241033309544)

and of course, you can get a sample of any size N For N=5: sum.getSample(5)

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