使用scipy.stats计算分布的平均值和标准差

3

我想要计算对数正态分布的均值和标准差,其中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,但是使用 lognorm.std(sigma,scale = np.exp(mu))仍然会得到 853.5698320847896 。


1
检查一下标准差的计算。当我修复你的代码时,我得到了500和另外一个数(165831.240)。 - cs95
1个回答

3

scipy.stats.lognorm 对数正态分布的参数化方式略有不同,为了与其他连续分布保持一致。第一个参数是形状参数,即您的 sigma。然后是 locscale 参数,它们允许对分布进行平移和缩放。在这里,您需要 loc=0.0scale=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

我手算得到的值与计算机结果一致。

为了确保这些参数选择确实给出了预期的对数正态分布,我创建了一个包含一百万个随机变量的样本,并查看了该样本的对数均值和标准差。如预期那样,这些值大致类似于你最初设定的musigma

>>> 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,您也可以通过解析的方式计算得到所需均值和标准差下的sigmamu。这样可以更快地获得更准确的结果:

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

谢谢你的帮助。我已经检查过了,发现有一些计算错误。现在我可以得到平均值=500,但仍然无法得到标准差=600。我已经将代码包含在编辑后的帖子中。 - VincentN
你的方差公式有误。请使用(exp(sigma**2) - 1)代替exp(sigma**2 - 1) - Mark Dickinson
它们之间有什么区别?因为它们仍然会产生相同的结果。 - VincentN
@VincentN:对我来说,使用exp(sigma**2) - 1可以得到正确的结果,而使用你原来的exp(sigma**2 - 1)则会得到错误的结果。至于“有什么区别”:exp(x - 1)exp(x) - 1是不同的函数,就像(x - 1)**2x**2 - 1是不同的函数一样。你为什么会期望它们是相同的呢? - Mark Dickinson
哦,我不知道昨晚为什么会忽略括号。现在它给出了正确的答案。谢谢你的帮助。 - VincentN

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接