计算截断对数正态分布的均值

Question

我正在尝试计算截断对数正态分布的平均值。我有一个随机变量x，它具有标准a的对数正态分布。

我想计算x < y时x的平均值

注意-如果x是正态分布的，可以使用this库来计算：

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)

我想转换这段代码，假设分布是对数正态分布，而不是正态分布。

Answer 1

也许有更好的方法可以做到这一点，但我最终还是回到了在截断结果之间的范围内对lognormal pdf乘以x进行积分来解决这个问题。

下面是一个Python示例-忽略我指定的未截断对数正态分布均值和标准差的笨拙方式，这只是我工作的一个特点。

它应该在任何截断(x1 =下限，x2 =上限)之间工作，包括0到无穷大(使用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=true&i=Integrate%5Bexp%5C%2840%29-Divide%5BPower%5B%5C%2840%29ln%5C%2840%29x%5C%2841%29-u%5C%2841%29%2C2%5D%2C%5C%2840%292*Power%5Bs%2C2%5D%5C%2841%29%5D%5C%2841%29%2Cx%5D
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