Блог

Я пытаюсь вычислить среднее значение усеченного логарифмического нормального распределения. У меня есть случайная величина x , которая имеет логарифмически нормальное распределение с std a .

Я хотел бы рассчитать среднее x значение, когда x < y

Примечание — Если x он был нормально распространен, его можно рассчитать с помощью этой библиотеки:

 from scipy.stats import truncnorm
my_mean = 100
my_std = 20
myclip_a = 0
myclip_b = 95
a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std
new_mean = truncnorm.mean(a, b, my_mean, my_std)
 

Я хотел бы преобразовать этот код в предположении, что распределение является логарифмически нормальным и ненормальным.

Вполне могут быть и более элегантные способы сделать это, но в итоге я вернулся к интеграции логнормального pdf-файла, умноженного на x в диапазоне между усеченными результатами, чтобы решить эту проблему раньше.

Ниже приведен пример Python — не обращайте внимания на неуклюжий способ, которым я указал среднее значение усеченного логнормального распределения и стандартное отклонение, это просто особенность моей работы.

Он должен работать между любыми усечениями (x1 = нижний предел, x2 = верхний предел), включая от нуля до бесконечности (с использованием np.inf)

 import math
from scipy.special import erf
import numpy as np

P10 = 50                         # Untruncated P10 (ie 10% outcomes higher than this)
P90 = 10                         # Untruncated P90  (ie 90% outcomes higher than this)
u = (np.log(P90) np.log(P10))/2  # Untruncated Mean of the log transformed distribution
s = np.log(P10/P90)/2.562        # Standard Deviation

# Returns integral of the lognormal pdf multiplied by the lognormal outcomes (x)
# Between lower (x1) and upper (x2) truncations
# pdf and cdf equations from https://en.wikipedia.org/wiki/Log-normal_distribution
# Integral evaluated with;
# https://www.wolframalpha.com/input/?i2d=trueamp;i=Integrate[exp(40)-Divide[Power[(40)ln(40)x(41)-u(41),2],(40)2*Power[s,2](41)](41),x]
def ln_trunc_mean(u, s, x1, x2):
    if x2 != np.inf:
        upper = erf((s**2 u-np.log(x2))/(np.sqrt(2)*s))
        upper_cum_prob = 0.5*(1 erf((np.log(x2)-u)/(s*np.sqrt(2)))) # CDF
    else:
        upper = -1
        upper_cum_prob = 1

    if x1 != 0:
        lower = erf((s**2 u-np.log(x1))/(np.sqrt(2)*s))
        lower_cum_prob = 0.5*(1 erf((np.log(x1)-u)/(s*np.sqrt(2))))
    else:
        lower = 1
        lower_cum_prob = 0

    integrand = -0.5*np.exp(s**2/2 u)*(upper-lower) # Integral of PDF.x.dx

    return integrand / (upper_cum_prob - lower_cum_prob)
 

Затем вы можете оценить, например, усеченное среднее значение, а также среднее значение с верхним и нижним 1 процентильным отсечением следующим образом

 # Untruncated Mean
print(ln_trunc_mean(u, s, 0, np.inf))
 

27.238164532490508

 # Truncated mean between 5.2 and 96.4
print(ln_trunc_mean(u, s, 5.2, 96.4))
 

26.5089880192863