Ottieni il valore della variabile casuale data la probabilità cumulativa (Python)

Aug 16 2020

Ecco una rapida informazione di base. Sto cercando di ottenere un CDF combinato per la combinazione lineare di due variabili casuali lognormali utilizzando l'approccio Monte-Carlo e quindi invertirlo per eseguire il campionamento. Ecco il codice Python per fare lo stesso:

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

Domande:

Esiste una funzione diretta in Scipy per eseguire questa operazione?

Nell'ultima riga del codice, sto prendendo il valore medio, c'è un modo per ottenere valori più accurati per interpolazione, ecc.? In tal caso, come lo implemento in Python

Risposte

3 SeverinPappadeux Aug 16 2020 at 23:36

Bene, c'è un caso ben noto in cui si sommano due RV X + Y, si conosce PDF X (x), PDF Y (y) e si desidera conoscere PDF X + Y (z). È possibile utilizzare un approccio simile qui, calcolare PDF e creare CDF = d PDF (z) / dz

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

dove Sdenota integrazione.

Potresti scriverlo direttamente per CDF

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

Puoi calcolare questo integrale:

  1. Analiticamente

  2. Numericamente, utilizzando SciPy

  3. Trasforma di Fourier in avanti e indietro, simile a Convolution

  4. Naturalmente, l'integrazione di Monte Carlo è sempre un'opzione

AGGIORNARE

Ecco il codice più semplice per iniziare

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

Se vuoi ottenere un campione dalla somma di 2 distribuzioni lognormali, non hai bisogno di uno schema 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 somma di x1 e x2 è essa stessa una distribuzione

sum = x1+x2

puoi accedere alla sua media sum.getMean()[0](= 1,5379) o alla sua deviazione standard sum.getStandardDeviation()[0](= 0,42689241033309544)

e, naturalmente, puoi ottenere un campione di qualsiasi dimensione N Per N = 5: sum.getSample(5)

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