python - 从python中的对数正态分布生成随机数

Question

我需要在 Python 中根据对数正态分布生成伪随机数。问题是我从对数正态 分布的众数 和标准差开始。我没有对数正态分布的均值或中位数，也没有基础正态分布的任何参数。

numpy.random.lognormal取基础正态分布的均值和标准差。我试图根据我拥有的参数来计算这些，但最终得到了一个四次函数。它有一个解决方案，但我希望有一个更直接的方法来做到这一点。

scipy.stats.lognorm接受我不明白的参数。我的母语不是英语，文档没有意义。

你能帮帮我吗？

Answer 1

您拥有对数正态分布的众数和标准差。要使用 scipy 的 lognorm 的 rvs() 方法，您必须根据形状参数 s 对分布进行参数化，这是标准的基础正态分布的偏差 sigma，以及 scale，即 exp(mu)，其中 mu 是基础分布的均值。

您指出进行此重新参数化需要求解四次多项式。为此，我们可以使用 numpy.poly1d 类。该类的实例具有 roots 属性。

一点代数表明 exp(sigma**2) 是多项式的唯一正实根

x**4 - x**3 - (stddev/mode)**2 = 0

其中 stddev 和 mode 是对数正态分布的给定标准差和模式，对于该解决方案，比例 (即 exp(mu)) 是

scale = mode*x

这是一个将模式和标准差转换为形状和尺度的函数:

def lognorm_params(mode, stddev):
    """
    Given the mode and std. dev. of the log-normal distribution, this function
    returns the shape and scale parameters for scipy's parameterization of the
    distribution.
    """
    p = np.poly1d([1, -1, 0, 0, -(stddev/mode)**2])
    r = p.roots
    sol = r[(r.imag == 0) & (r.real > 0)].real
    shape = np.sqrt(np.log(sol))
    scale = mode * sol
    return shape, scale

例如，

In [155]: mode = 123

In [156]: stddev = 99

In [157]: sigma, scale = lognorm_params(mode, stddev)

使用计算的参数生成样本:

In [158]: from scipy.stats import lognorm

In [159]: sample = lognorm.rvs(sigma, 0, scale, size=1000000)

这是样本的标准差:

In [160]: np.std(sample)
Out[160]: 99.12048952171304

下面是一些 matplotlib 代码，用于绘制样本的直方图，并在抽取样本的分布模式处绘制一条垂直线:

In [176]: tmp = plt.hist(sample, normed=True, bins=1000, alpha=0.6, color='c', ec='c')

In [177]: plt.xlim(0, 600)
Out[177]: (0, 600)

In [178]: plt.axvline(mode)
Out[178]: <matplotlib.lines.Line2D at 0x12c5a12e8>

直方图:

---

如果您想使用 numpy.random.lognormal() 而不是 scipy.stats.lognorm.rvs() 生成样本，您可以这样做:

In [200]: sigma, scale = lognorm_params(mode, stddev)

In [201]: mu = np.log(scale)

In [202]: sample = np.random.lognormal(mu, sigma, size=1000000)

In [203]: np.std(sample)
Out[203]: 99.078297384090902

---

我没有研究过 poly1d 的 roots 算法有多稳健，所以一定要测试各种可能的输入值。或者，您可以使用 scipy 中的求解器来求解上述 x 的多项式。您可以使用以下方式绑定(bind)解决方案:

max(sqrt(stddev/mode), 1) <= x <= sqrt(stddev/mode) + 1