时间戳

首页 > 解决方案 > 计算截断对数正态分布的平均值

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

问题描述

我正在尝试计算截断对数正态分布的平均值。我有一个随机变量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)

我想在假设分布是对数正态而不是正态的情况下转换此代码。

标签： mathscipystatistics

解决方案

可能有更优雅的方法可以做到这一点，但我最终恢复到在截断结果之间的范围内整合对数正态 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=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

推荐阅读