Ermitteln Sie den Wert der Zufallsvariablen unter Berücksichtigung der kumulativen Wahrscheinlichkeit (Python).

Aug 16 2020

Hier finden Sie eine kurze Hintergrundinformation. Ich versuche, eine kombinierte CDF für die lineare Kombination zweier lognormaler Zufallsvariablen unter Verwendung des Monte-Carlo-Ansatzes zu erhalten und sie dann zu invertieren, um eine Stichprobe zu erstellen. Hier ist der Python-Code, um dasselbe zu tun:

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

Fragen:

Gibt es in scipy eine direkte Funktion, um diesen Vorgang auszuführen?

In der letzten Zeile des Codes nehme ich den Mittelwert. Gibt es eine Möglichkeit, durch Interpolation usw. genauere Werte zu erhalten? Wenn ja, wie implementiere ich es in Python?

Antworten

3 SeverinPappadeux Aug 16 2020 at 23:36

Nun, es gibt einen bekannten Fall, wenn Sie zwei Wohnmobile X + Y summieren, PDF X (x), PDF Y (y) kennen und PDF X + Y (z) kennen wollen. Sie können hier einen ähnlichen Ansatz verwenden, PDF berechnen und CDF = d PDF (z) / dz erstellen

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

wo Sbezeichnet Integration.

Sie können es direkt für CDF schreiben

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

Sie könnten dieses Integral berechnen:

  1. Analytisch

  2. Numerisch mit SciPy

  3. Führen Sie eine Fourier-Transformation vorwärts und rückwärts durch, ähnlich wie bei der Faltung

  4. Natürlich ist die Monte-Carlo-Integration immer eine Option

AKTUALISIEREN

Hier ist der einfachste Code, mit dem Sie loslegen können

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

JeanA. Oct 22 2020 at 00:55

Wenn Sie eine Stichprobe aus der Summe von 2 logarithmischen Normalverteilungen erhalten möchten, benötigen Sie kein Monte-Carlo-Schema.

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

Die Summe von x1 und x2 ist selbst eine Verteilung

sum = x1+x2

Sie können auf den Mittelwert sum.getMean()[0](= 1,5379) oder die Standardabweichung sum.getStandardDeviation()[0](= 0,42689241033309544) zugreifen.

und natürlich können Sie eine Stichprobe jeder Größe N erhalten. Für N = 5: sum.getSample(5)

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