Obtenga el valor de la variable aleatoria dada la probabilidad acumulada (Python)

Aug 16 2020

Aquí hay una rápida información de antecedentes. Estoy tratando de obtener un CDF combinado para la combinación lineal de dos variables aleatorias logarítmicas normales utilizando el enfoque de Monte-Carlo y luego invertirlo para hacer un muestreo. Aquí está el código de Python para hacer lo mismo:

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

Preguntas:

¿Existe una función directa en scipy para realizar esta operación?

En la última línea del código, estoy tomando el valor medio, ¿hay alguna manera de obtener valores más precisos por interpolación, etc.? Si es así, ¿cómo lo implemento en Python?

Respuestas

3 SeverinPappadeux Aug 16 2020 at 23:36

Bueno, hay un caso bien conocido en el que sumas dos RVs X + Y, conoces PDF X (x), PDF Y (y) y quieres saber PDF X + Y (z). Puede usar un enfoque similar aquí, calcular PDF y hacer CDF = d PDF (z) / dz

PDF aX + bY (z) = S dy PDF Y (y) PDF X ((z-by) / a) / | a |

donde Sdenota integración.

Podrías escribirlo directamente para CDF

CDF aX + bY (z) = S dy PDF Y (y) CDF X ((z-by) / a)

Podrías calcular esta integral:

  1. Analíticamente

  2. Numéricamente, usando SciPy

  3. ¿Transformación de Fourier hacia adelante y hacia atrás, similar a la convolución?

  4. Por supuesto, la integración de Monte Carlo es siempre una opción

ACTUALIZAR

Aquí está el código más simple para comenzar

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

Grafico

JeanA. Oct 22 2020 at 00:55

Si desea obtener una muestra de la suma de 2 distribuciones logarítmicas normales, no necesita un esquema de 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]))

la suma de x1 y x2 es en sí misma una distribución

sum = x1+x2

puede acceder a su media sum.getMean()[0](= 1,5379) o su desviación estándar sum.getStandardDeviation()[0](= 0,42689241033309544)

y, por supuesto, puede obtener una muestra de cualquier tamaño N Para N = 5: sum.getSample(5)

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