如何在Python中使用Mu和Sigma获得对数正态分布?

34

我一直在尝试使用Scipy获取lognormal分布的结果。 我已经拥有Mu和Sigma,所以我不需要做任何其他准备工作。如果我需要更具体(而我正试图用有限的统计知识变得更具体),我会说我正在寻找累积函数(Scipy下的cdf)。问题是,我无法弄清楚如何仅使用0-1范围内的平均值和标准偏差来完成这项任务(即返回的答案应该是0-1的某些东西)。我也不确定应该使用dist 中的哪种方法来获得答案。我尝试阅读文档和查看SO,但相关问题(例如thisthis)似乎没有提供我要寻找的答案。

这里是我正在使用的代码示例。谢谢。

from scipy.stats import lognorm
stddev = 0.859455801705594
mean = 0.418749176686875
total = 37
dist = lognorm.cdf(total,mean,stddev)

更新:

经过一些工作和一点研究,我有了一些进展。但是我仍然得到了错误的答案。新代码如下。根据 R 和 Excel,结果应该是0.7434,但显然不是这样的。我是否遗漏了逻辑上的缺陷?

dist = lognorm([1.744],loc=2.0785)
dist.cdf(25)  # yields=0.96374596, expected=0.7434

更新2: 工作日志正态实现,可得到正确的0.7434结果。

def lognorm(self,x,mu=0,sigma=1):
   a = (math.log(x) - mu)/math.sqrt(2*sigma**2)
   p = 0.5 + 0.5*math.erf(a)
   return p
lognorm(25,1.744,2.0785)
> 0.7434

2
你能解释一下你对“分布结果”的理解吗? - joaquin
@joaquin 我添加了一个代码示例,展示了我拥有的内容以及我期望它产生的结果。 - Eric Lubow
@EricLubow:我认为你可能误解了这种情况下mean和stddev的含义。对于对数正态分布,它们是变量对数的平均值和标准差。如果一个变量是对数正态分布的,那么它意味着该变量的对数是正态分布的。 - talonmies
@talonmies 我了解使用平均数和标准差意味着使用变量的对数的平均数和标准差。我手写了Python的lognorm函数(如上所示)并且能够得到正确的答案。这就是让我相信在SciPy实现中可能存在差异的原因,因为我在R和Excel中得到了正确的答案。如果我的实现有误,我一定想知道。 - Eric Lubow
@EricLubow:你基本上重新实现了scipy内部使用的内容。Lucas的答案似乎是正确的。请参见我的答案以获取用法示例。 - serv-inc
对于新读者发现这个问题 - 我发现这个问题的被接受的答案比下面的答案更有用。 - charelf
7个回答

46

我知道这有点晚了(快一年了!),但我一直在研究scipy.stats中的lognorm函数。很多人似乎对输入参数感到困惑,所以我希望能帮助这些人。上面的示例几乎是正确的,但我发现将均值设置为“loc”参数有点奇怪——这表明cdf或pdf在值大于均值时才开始“起飞”。此外,均值和标准差参数应分别以exp(Ln(mean))和Ln(StdDev)的形式给出。

简而言之,参数为(x,形状,loc,scale),其中参数定义如下:

loc-无等效项,这将从您的数据中减去,使0成为数据范围的infimum。

scale- exp μ,其中μ是变量的对数的平均值。(当拟合时,通常使用数据对数的样本平均值。)

形状-变量的对数的标准偏差。

我经历了与大多数人一样的挫败感,所以我分享了我的解决方案。只要小心,因为没有资源汇编,这些解释就不太清楚。

更多信息,请参见以下有用的来源:

以下是示例,取自@serv-inc发布在此页面here的答案:

import math
from scipy import stats

# standard deviation of normal distribution
sigma = 0.859455801705594
# mean of normal distribution
mu = 0.418749176686875
# hopefully, total is the value where you need the cdf
total = 37

frozen_lognorm = stats.lognorm(s=sigma, scale=math.exp(mu))
frozen_lognorm.cdf(total) # use whatever function and value you need here

3
如果我理解正确的话:在数学符号中,如果X服从N(mu,sigma)分布,则Y=exp(X)服从LogN(mu,sigma)分布。要在Scipy中获取X,可以使用norm(mu,sigma),但要获取Y,则需要使用lognorm(sigma,0,exp(mu))。这有些棘手... - Elmar Zander
5
顺便说一句,我觉得你的帖子很有帮助,但是scipy的帮助文档就不太行了。对于每个分布,你都需要尝试一下参数的含义(例如对于均匀分布U(a,b),其中[a,b]是区间,你需要使用uniform(loc=a, scale=b-a)来生成,这里的loc不是均值,scale也不是标准差...)。 - Elmar Zander
4
дҪ еҸҜд»ҘдҪҝз”Ёlognorm(s=sigma, scale=math.exp(mu)пјҢиҜҰжғ…иҜ·и§Ғhttps://dev59.com/questions/Cmox5IYBdhLWcg3w8I1J#36714419 - serv-inc
1
@Lucas 晚来总比不来好;) 非常感谢。非常有用的总结。 - Aleksander Lidtke

24

听起来你想从已知参数实例化一个“冻结”的分布。在你的例子中,你可以这样做:

from scipy.stats import lognorm
stddev = 0.859455801705594
mean = 0.418749176686875
dist=lognorm([stddev],loc=mean)

使用这个方法可以得到一个带有指定均值和标准差的对数正态分布对象。然后,您可以像下面这样获取pdf或cdf:

import numpy as np
import pylab as pl
x=np.linspace(0,6,200)
pl.plot(x,dist.pdf(x))
pl.plot(x,dist.cdf(x))

对数正态分布的累积分布函数和概率密度函数

这是你想要的吗?


1
笔误:应该是"np.linspace"而不是"np.inspace"。 - Max Li
应该不是“dist=lognorm([stddev**2],loc=mean)”,而是方差作为参数,而不是标准差吗?我在scipy文档中没有找到参数规范,你知道吗? - Max Li
我更新了我的问题,展示了我目前遇到的问题。我认为问题可能在于SciPy对于对数正态分布函数(cdf)的实现与R或Excel不同。当然,也有可能是我使用了错误的方法。 - Eric Lubow
2
根据下面Lucas的回答,这是错误的,对吗?平均值不应该在分布的最左边,而应该在峰值右侧,对吗? - Alex S
@AlexS:Lucas 说得对。有关代码示例,请参见 https://dev59.com/Cmox5IYBdhLWcg3w8I1J#36714419 - serv-inc
显示剩余2条评论

14
from math import exp
from scipy import stats

def lognorm_cdf(x, mu, sigma):
    shape  = sigma
    loc    = 0
    scale  = exp(mu)
    return stats.lognorm.cdf(x, shape, loc, scale)

x      = 25
mu     = 2.0785
sigma  = 1.744
p      = lognorm_cdf(x, mu, sigma)  #yields the expected 0.74341

与Excel和R类似,上面的lognorm_cdf函数使用musigma参数对对数正态分布进行累积分布函数(CDF)参数化。

虽然SciPy使用shapelocscale参数来表征其概率分布,但是对于对数正态分布,我发现将这些参数视为变量级别而不是分布级别会稍微容易一些。这就是我的意思...

一个对数正态变量X与一个正态变量Z的关系如下:

X = exp(mu + sigma * Z)              #Equation 1

这与以下内容相同:

X = exp(mu) * exp(Z)**sigma          #Equation 2

这个可以巧妙地重写如下:

X = exp(mu) * exp(Z-Z0)**sigma       #Equation 3

其中Z0 = 0。这个方程的形式为:

f(x) = a * ( (x-x0) ** b )           #Equation 4

如果你能在脑海中形象地想像方程,那么方程4中的比例、形状和位置参数分别是:abx0。 这意味着在方程3中,比例、形状和位置参数分别为:exp(mu)sigma和零。

如果你不能非常清楚地想象出它,让我们将方程2重新写成一个函数:

f(Z) = exp(mu) * exp(Z)**sigma      #(same as Equation 2)

接下来我们看一下musigmaf(Z)的影响。下面的图像保持sigma不变,改变mu。您应该可以看到,mu会垂直缩放f(Z)。然而,它并不是线性变换;将mu从0变为1所产生的效果比将mu从1变为2产生的效果要小。根据方程2,我们可以看出exp(mu)实际上是线性缩放因子。因此SciPy的"scale"就是exp(mu)

effects_of_mu

下一个图像保持mu不变,改变sigma。您应该可以看到f(Z)的形状发生了变化。也就是说,当Z=0时,f(Z)具有恒定值,而sigma影响f(Z)沿水平轴迅速向曲线弯曲。因此SciPy的"shape"就是sigma

effects_of_sigma


能否解释一下为什么这是问题的答案? - MeanGreen
我发现这与Excel函数LOGNORM.DIST(x,Mu,Sigma,TRUE)一一对应。 - asdag8

4

@lucas' answer已经很好地解释了使用方法。作为代码示例,你可以使用

import math
from scipy import stats

# standard deviation of normal distribution
sigma = 0.859455801705594
# mean of normal distribution
mu = 0.418749176686875
# hopefully, total is the value where you need the cdf
total = 37

frozen_lognorm = stats.lognorm(s=sigma, scale=math.exp(mu))
frozen_lognorm.cdf(total) # use whatever function and value you need here

3

已知对数正态分布的平均值和标准差

如果有人需要获取 scipy.stats.lognorm 分布,并且已知对数正态分布的平均值 mu 和标准差 sigma,那么这里有一个解决方案。在这种情况下,我们需要从已知的 musigma 计算出 stats.lognorm 参数,方法如下:

import numpy as np
from scipy import stats

mu = 10
sigma = 3

a = 1 + (sigma / mu) ** 2
s = np.sqrt(np.log(a))
scale = mu / np.sqrt(a)

这是通过查看stats.lognorm.stats方法中方差和均值计算的实现方式,并基本上反向操作(解决输入)获得的。

然后,我们可以初始化冻结的分布实例。

distr = stats.lognorm(s, 0, scale)

# generate some randomvals
randomvals = distr.rvs(1_000_000)
# calculate mean and variance using the dedicated method
mu_stats, var_stats = distr.stats("mv")

比较输入、随机值和分布解析解的均值和标准差,可使用 distr.stats 中的函数。

print(f"""
                 Mean    Std
----------------------------
Input:         {mu:6.2f} {sigma:6.2f}
Randomvals:    {randomvals.mean():6.2f} {randomvals.std():6.2f}
lognorm.stats: {mu_stats:6.2f} {np.sqrt(var_stats):6.2f}
""")

                 Mean    Std
----------------------------
Input:          10.00   3.00
Randomvals:     10.00   3.00
lognorm.stats:  10.00   3.00

使用stats.lognorm绘制概率密度函数图和随机值的直方图:

import holoviews as hv
hv.extension('bokeh')

x = np.linspace(0, 30, 301)
counts, _ = np.histogram(randomvals, bins=x)
counts = counts / counts.sum() / (x[1] - x[0])

(hv.Histogram((counts, x)) 
* hv.Curve((x, distr.pdf(x))).opts(color="r").opts(width=900))

enter image description here


3
即使有点晚了,但如果对其他人有帮助:我发现 Excel 的 TEXT() 函数可以将数字格式化为文本。
LOGNORM.DIST(x,Ln(mean),standard_dev,TRUE)

提供与Python相同的结果。
from scipy.stats import lognorm
lognorm.cdf(x,sigma,0,mean)

同样地,Excel 的
LOGNORM.DIST(x,Ln(mean),standard_dev,FALSE)

似乎等同于Python的
from scipy.stats import lognorm
lognorm.pdf(x,sigma,0,mean).

对于第一个案例,在 x=2039.9337873,mean=7.6901,std_dev=0.6772 的情况下,它们没有给我返回相同的结果。 - Deniz Ozger
啊,我忘记在我的Excel公式中添加Ln(mean)了。已在答案中进行更正。 - Docuemada

1
如果您只是想要一个行为类似于R中的lnorm函数的功能。那么,请放下暴怒,使用numpy的numpy.random.lognormal

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