Obtenir la valeur de la variable aléatoire compte tenu de la probabilité cumulée (Python)

Aug 16 2020

Voici une brève information de base. J'essaie d'obtenir un CDF combiné pour la combinaison linéaire de deux variables aléatoires log-normales en utilisant l'approche de Monte-Carlo, puis de l'inverser pour faire de l'échantillonnage. Voici le code Python pour faire de même:

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

Des questions:

Existe-t-il une fonction directe dans scipy pour effectuer cette opération?

Dans la dernière ligne du code, je prends la valeur moyenne, y a-t-il un moyen d'obtenir des valeurs plus précises par interpolation, etc.? Si oui, comment l'implémenter en Python

Réponses

3 SeverinPappadeux Aug 16 2020 at 23:36

Eh bien, il existe un cas bien connu lorsque vous additionnez deux RV X + Y, connaissez PDF X (x), PDF Y (y) et voulez connaître PDF X + Y (z). Vous pouvez utiliser une approche similaire ici, calculer un PDF et créer CDF = d PDF (z) / dz

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

Sdénote l'intégration.

Vous pouvez l'écrire directement pour CDF

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

Vous pouvez calculer cette intégrale:

  1. Analytiquement

  2. Numériquement, en utilisant SciPy

  3. Faire une transformation de Fourier en avant et en arrière, similaire à la convolution

  4. Bien sûr, l'intégration Monte Carlo est toujours une option

METTRE À JOUR

Voici le code le plus simple pour vous aider

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

Graphique

JeanA. Oct 22 2020 at 00:55

Si vous souhaitez obtenir un échantillon de la somme de 2 distributions log-normales, vous n'avez pas besoin d'un schéma 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 somme de x1 et x2 est elle-même une distribution

sum = x1+x2

vous pouvez accéder à sa moyenne sum.getMean()[0](= 1,5379) ou à son écart type sum.getStandardDeviation()[0](= 0,42689241033309544)

et bien sûr, vous pouvez obtenir un échantillon de n'importe quelle taille N Pour N = 5: sum.getSample(5)

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