使用 scipy.stats 分布均值和标准差

Question

我试图获得对数正态分布的均值和标准差，其中 mu=0.4104857306 和 sigma=3.4070874277012617，我期望均值 = 500 和标准差 = 600。我不确定我做错了什么。代码如下：

import scipy.stats as stats
import numpy as np
a = 3.4070874277012617
b = 0.4104857306
c = stats.lognorm.mean(a,b)
d = stats.lognorm.var(a,b)
e = np.sqrt(d)
print("Mean =",c)
print("std =",e)

输出在这里：

Mean = 332.07447304207903
sd = 110000.50047821256

提前谢谢你。

编辑：

感谢您的帮助。我检查并发现有一些计算错误。我现在可以得到平均值 = 500，但仍然无法得到标准 = 600。这是我使用的代码：

import numpy as np
import math
from scipy import exp
from scipy.optimize import fsolve

def f(z):
    mean = 500
    std = 600
    sigma = z[0]
    mu = z[1]
    f = np.zeros(2)
    f[0] = exp(mu + (sigma**2) / 2) - mean
    f[1] = exp(2*mu + sigma**2) * exp(sigma**2 - 1) - std**2
    return f
z = fsolve (f,[1.1681794012855686,5.5322865416282365])
print("sigma =",z[0])
print("mu =",z[1])
print(f(z))

sigma = 1.1681794012855686
mu = 5.5322865416282365

我试过用我的计算器检查结果，我可以按要求得到std=600，我仍然得到853.5698320847896和lognorm.std(sigma, scale=np.exp(mu))。

Answer 1

scipy.stats.lognorm 对数正态分布以一种稍微不寻常的方式进行参数化，以便与其他连续分布保持一致。第一个参数是形状参数，即您的sigma。紧随其后的是loc 和scale 参数，它们允许对分布进行移动和缩放。在这里你想要loc=0.0 和scale=exp(mu)。因此，要计算平均值，您需要执行以下操作：

>>> import numpy as np
>>> from scipy.stats import lognorm
>>> mu = 0.4104857306
>>> sigma = 3.4070874277012617
>>> lognorm.mean(sigma, 0.0, np.exp(mu))
500.0000010889041

或者更清楚一点：按名称传递scale 参数，并将loc 参数保留为默认的0.0：

>>> lognorm.mean(sigma, scale=np.exp(mu))
500.0000010889041

正如@coldspeed 在他的评论中所说，您对标准差的预期值看起来不正确。我明白了：

>>> lognorm.std(sigma, scale=np.exp(mu))
165831.2402402415

我手动计算得到相同的值。

为了仔细检查这些参数选择是否确实给出了预期的对数正态分布，我创建了一个包含一百万个偏差的样本，并查看了该样本对数的均值和标准偏差。正如预期的那样，这些返回的值与您原来的 mu 和 sigma 大致相似：

>>> samples = lognorm.rvs(sigma, scale=np.exp(mu), size=10**6)
>>> np.log(samples).mean()  # should be close to mu
0.4134644116056518
>>> np.log(samples).std(ddof=1)  # should be close to sigma
3.4050012251732285

---

响应编辑：您的对数正态方差公式略有错误 - 您需要将 exp(sigma**2 - 1) 替换为 (exp(sigma**2) - 1)。如果你这样做，并重新运行 fsolve 计算，你会得到：

sigma = 0.9444564779275075
mu = 5.768609079062494

使用这些值，您应该得到预期的均值和标准差：

>>> from scipy.stats import lognorm
>>> import numpy as np
>>> sigma = 0.9444564779275075
>>> mu = 5.768609079062494
>>> lognorm.mean(sigma, scale=np.exp(mu))
499.9999999949592
>>> lognorm.std(sigma, scale=np.exp(mu))
599.9999996859631

除了使用fsolve，您还可以解析求解sigma 和mu，给定所需的均值和标准差。这可以更快地为您提供更准确的结果：

>>> mean = 500.0
>>> var = 600.0**2
>>> sigma = np.sqrt(np.log1p(var/mean**2))
>>> mu = np.log(mean) - 0.5*sigma*sigma
>>> mu, sigma
(5.768609078769636, 0.9444564782482624)
>>> lognorm.mean(sigma, scale=np.exp(mu))
499.99999999999966
>>> lognorm.std(sigma, scale=np.exp(mu))
599.9999999999995