Scipy，对数正态分布-参数

Question

Scipy，对数正态分布-参数

pythonstatisticsscipy

39

我想使用Python的`scipy.stats.lognormal.fit`函数将我的数据拟合成对数正态分布。根据手册，`fit`返回形状，位置，比例参数。但是，对数正态分布通常只需要两个参数：均值和标准差。

如何解释scipy `fit`函数的结果？怎样得到均值和标准差？

- Jakub M.

5个回答

15

我刚刚花了一些时间研究这个问题，并想在这里记录一下：如果你想从lognorm.fit的三个返回值 ((shape, loc, scale)) 中获取概率密度 (在点x)，你需要使用以下公式：

x = 1 / (shape*((x-loc)/scale)*sqrt(2*pi)) * exp(-1/2*(log((x-loc)/scale)/shape)**2) / scale

因此，给定方程式为(loc为µ，shape为σ，scale为α)：

$x = \frac{1}{(x-\mu)\cdot\sqrt{2\pi\sigma^2}} \cdot e^{-\frac{log(\frac{x-\mu}{\alpha})^2}{2\sigma^2}}$

- Chronial

你的公式末尾为什么要除以(1/alpha)呢？ - bioslime

2

这不是我的公式，而是scipy的工作方式。我没有看到最后除以（1/α}，所以我会假设您说的是除以scale - 如果我误解了，请纠正我。我一段时间前看过这段代码，但如果我记得正确，那就是对每种分布都要做的事情。但请注意，在开始时还有另一个scale，它们互相抵消了（正如您在方程中看到的那样）。 - Chronial

这两个比例尺相互抵消。所以最好使用它： (np.exp((-1/2)*(np.log((x-loc)/scale)/shape)*2))/(shape(x-loc)np.sqrt(2np.pi)) - Gouz

2

我认为这会有所帮助。 我长时间以来一直在寻找同样的问题，并最终找到了解决我的问题的方法。在我的情况下，我试图使用scipy.stats.lognorm模块将一些数据适配到对数正态分布上。 然而，当我最终得到模型参数时，我无法找到一种方法使用y数据的平均值和标准差复制我的结果。

在下面的代码中，我从平均值和标准差参数中解释如何使用scipy.stats.norm模块生成正态分布的数据样本。使用这些数据，我适配了正常模型（norm_dist_fitted），并创建了一个使用从数据中提取的平均值和标准偏差（mu，sigma）的正常模型。

原始模型产生的数据，拟合和由（mu-sigma）对产生的数据在图表中进行了比较。

Fig1

在代码的下一部分，我使用普通数据生成对数正态分布的样本。为此，请注意对数正态样本将是原始样本的指数。因此，指数样本的均值和标准差将分别为（exp(mu) 和 exp(sigma)）。

我将生成的数据拟合到一个lognormal模型中（由于我的样本的对数（exp(x)) 是正态分布的，并且遵循对数正态模型的假设）。

要从原始数据（x）的均值和标准差生成对数正态模型，代码如下：

lognorm_dist = scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu))

然而，如果你的数据已经在指数空间（exp(x)）中，则必须使用：

muX = np.mean(np.log(x))
sigmaX = np.std(np.log(x))
scipy.stats.lognorm(s=sigmaX, loc=0, scale=muX)

import scipy
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

mu = 10 # Mean of sample !!! Make sure your data is positive for the lognormal example 
sigma = 1.5 # Standard deviation of sample
N = 2000 # Number of samples

norm_dist = scipy.stats.norm(loc=mu, scale=sigma) # Create Random Process
x = norm_dist.rvs(size=N) # Generate samples

# Fit normal
fitting_params = scipy.stats.norm.fit(x)
norm_dist_fitted = scipy.stats.norm(*fitting_params)
t = np.linspace(np.min(x), np.max(x), 100)

# Plot normals
f, ax = plt.subplots(1, sharex='col', figsize=(10, 5))
sns.distplot(x, ax=ax, norm_hist=True, kde=False, label='Data X~N(mu={0:.1f}, sigma={1:.1f})'.format(mu, sigma))
ax.plot(t, norm_dist_fitted.pdf(t), lw=2, color='r',
        label='Fitted Model X~N(mu={0:.1f}, sigma={1:.1f})'.format(norm_dist_fitted.mean(), norm_dist_fitted.std()))
ax.plot(t, norm_dist.pdf(t), lw=2, color='g', ls=':',
        label='Original Model X~N(mu={0:.1f}, sigma={1:.1f})'.format(norm_dist.mean(), norm_dist.std()))
ax.legend(loc='lower right')
plt.show()


# The lognormal model fits to a variable whose log is normal
# We create our variable whose log is normal 'exponenciating' the previous variable

x_exp = np.exp(x)
mu_exp = np.exp(mu)
sigma_exp = np.exp(sigma)

fitting_params_lognormal = scipy.stats.lognorm.fit(x_exp, floc=0, scale=mu_exp)
lognorm_dist_fitted = scipy.stats.lognorm(*fitting_params_lognormal)
t = np.linspace(np.min(x_exp), np.max(x_exp), 100)

# Here is the magic I was looking for a long long time
lognorm_dist = scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu))

# The trick is to understand these two things:
# 1. If the EXP of a variable is NORMAL with MU and STD -> EXP(X) ~ scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu))
# 2. If your variable (x) HAS THE FORM of a LOGNORMAL, the model will be scipy.stats.lognorm(s=sigmaX, loc=0, scale=muX)
# with:
#    - muX = np.mean(np.log(x))
#    - sigmaX = np.std(np.log(x))


# Plot lognormals
f, ax = plt.subplots(1, sharex='col', figsize=(10, 5))
sns.distplot(x_exp, ax=ax, norm_hist=True, kde=False,
             label='Data exp(X)~N(mu={0:.1f}, sigma={1:.1f})\n X~LogNorm(mu={0:.1f}, sigma={1:.1f})'.format(mu, sigma))
ax.plot(t, lognorm_dist_fitted.pdf(t), lw=2, color='r',
        label='Fitted Model X~LogNorm(mu={0:.1f}, sigma={1:.1f})'.format(lognorm_dist_fitted.mean(), lognorm_dist_fitted.std()))
ax.plot(t, lognorm_dist.pdf(t), lw=2, color='g', ls=':',
        label='Original Model X~LogNorm(mu={0:.1f}, sigma={1:.1f})'.format(lognorm_dist.mean(), lognorm_dist.std()))
ax.legend(loc='lower right')
plt.show()

- nenetto

请在回答中提供一些解释。 - Collin Barrett

@CollinM.Barrett 完成了！ :D 希望一切都好 :D - nenetto

0

首先，loc 不是分布的简单线性偏移，实际上，loc 有其自身的统计意义，它表示样本减去 loc 后会得到一个“标准化”的对数正态分布，其下限为零，这一点非常重要。

因此，当您指定“loc”或“floc”时，实际上您施加了一种非常强的假设，即您假定这些样本具有下限，并且下限“恰好”是“loc”的值。因此，scipy 使用不同的算法进行拟合，即：如果您提供 loc 信息，则 scipy 将采用最大似然方法计算拟合参数；如果没有，则将使用数值解算器。

此外，您可以查看代码：在 scipy 包 stats/_continuous_distns.py 第 3889 行。如下所示：

    def fit(self, data, *args, **kwds):
        floc = kwds.get('floc', None)
        if floc is None:
            # loc is not fixed.  Use the default fit method.
            return super(lognorm_gen, self).fit(data, *args, **kwds)

        f0 = (kwds.get('f0', None) or kwds.get('fs', None) or
              kwds.get('fix_s', None))
        fscale = kwds.get('fscale', None)

        if len(args) > 1:
            raise TypeError("Too many input arguments.")
        for name in ['f0', 'fs', 'fix_s', 'floc', 'fscale', 'loc', 'scale',
                     'optimizer']:
            kwds.pop(name, None)
        if kwds:
            raise TypeError("Unknown arguments: %s." % kwds)

        # Special case: loc is fixed.  Use the maximum likelihood formulas
        # instead of the numerical solver.

此外，来自R社区的人可能会想知道为什么python的输出与R不同。实际上，我不赞成将R用作“参考”，它只是一种软件，不同的软件有不同风味的算法。

例如，R的输出如下所示并非错误，Python或其他软件（如Fortran）使用完全不同的算法：

round(3.5)
[1] 4
round(2.5)
[1] 2

- Chao Sun

-1

有助于我的一点是将位置和比例作为参数化考虑。

与标准对数正态分布中使用x不同，您可以更改为x'=(x-位置)/比例

然后概率密度函数F(x')=(1/比例)F((x-位置)/比例))

更多信息请参见链接https://en.wikipedia.org/wiki/Location%E2%80%93scale_family

- Mateus Kinasz

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- ars · Accepted Answer

在scipy中，分布是以通用的方式编码的，涉及两个参数：位置和尺度。其中位置参数（loc）将分布向左或向右移动，而尺度参数则压缩或拉伸分布。对于双参数对数正态分布，"均值"和"标准差"分别对应于log(scale)和shape（您可以让loc=0）。以下演示了如何拟合对数正态分布以找到感兴趣的两个参数。

In [56]: import numpy as np

In [57]: from scipy import stats

In [58]: logsample = stats.norm.rvs(loc=10, scale=3, size=1000) # logsample ~ N(mu=10, sigma=3)

In [59]: sample = np.exp(logsample) # sample ~ lognormal(10, 3)

In [60]: shape, loc, scale = stats.lognorm.fit(sample, floc=0) # hold location to 0 while fitting

In [61]: shape, loc, scale
Out[61]: (2.9212650122639419, 0, 21318.029350592606)

In [62]: np.log(scale), shape  # mu, sigma
Out[62]: (9.9673084420467362, 2.9212650122639419)