正态分布与标准差不匹配

Question

作为研究的一部分，我测量了对数正态分布的平均值和标准差。给定一个基本的正态分布值，就可以分析性地预测这些量（如https://en.wikipedia.org/wiki/Log-normal_distribution）了。在

然而，从下面的情节可以看出，情况似乎并非如此。第一个曲线图显示对数正态数据相对于高斯分布的sigma的平均值，而第二个图显示对数正态数据相对于高斯sigma的sigma。显然，“计算”线与“分析”线的偏离非常明显。在

我把高斯分布的平均值与mu = -0.5*sigma**2联系起来，因为这确保了对数正态场的平均值应该是1。注意，这是由我工作的物理领域引起的：例如，如果设置mu=0.0，那么与解析值的偏差仍然会发生。在

通过复制并粘贴问题底部的代码，可以复制下面的图。任何关于可能导致这种情况的建议都将不胜感激！在

对数正态与高斯西格玛的平均值：

对数正态的西格玛与高斯的西格玛：

注意，为了生成上面的图，我使用了N=10000，但是为了速度，在下面的代码中添加了{}。在

import numpy as np
import matplotlib.pyplot as plt

mean_calc = []
sigma_calc = []
mean_analytic = []
sigma_analytic = []
ss = np.linspace(1.0,10.0,46)
N = 1000

for s in ss:
  mu = -0.5*s*s
  ln = np.random.lognormal(mean=mu, sigma=s, size=(N,N))
  mean_calc += [np.average(ln)]
  sigma_calc += [np.std(ln)]
  mean_analytic += [np.exp(mu+0.5*s*s)]
  sigma_analytic += [np.sqrt((np.exp(s**2)-1)*(np.exp(2*mu + s*s)))]

plt.loglog(ss,mean_calc,label='calculated')
plt.loglog(ss,mean_analytic,label='analytic')
plt.legend();plt.grid()
plt.xlabel(r'$\sigma_G$')
plt.ylabel(r'$\mu_{LN}$')
plt.show()

plt.loglog(ss,sigma_calc,label='calculated')
plt.loglog(ss,sigma_analytic,label='analytic')
plt.legend();plt.grid()
plt.xlabel(r'$\sigma_G$')
plt.ylabel(r'$\sigma_{LN}$')
plt.show()

Answer 1

TL；DR

对数正态分布为正偏态和重尾分布。当对从高度偏态分布中提取的样本执行浮点算术运算（如sum、mean或std）时，采样向量包含的值在几个数量级（几十年）内存在差异。这使得计算不准确。在

问题来自这两条线：

mean_calc += [np.average(ln)]
sigma_calc += [np.std(ln)]

因为ln包含非常低和非常高的值，其数量级远远高于浮点精度。在

使用以下谓词可以很容易地检测到问题，以警告用户其计算错误：

^{pr2}$

这在严格正实中显然是错误的，但在使用有限精度算法时必须加以考虑。在

修改你的MCVE

如果我们稍微修改一下您的MCVE为：

from scipy import stats

for s in ss:
    mu = -0.5*s*s
    ln = stats.lognorm(s, scale=np.exp(mu)).rvs(N*N)
    f = stats.lognorm.fit(ln, floc=0)
    mean_calc += [f[2]*np.exp(0.5*s*s)]
    sigma_calc += [np.sqrt((np.exp(f[0]**2)-1)*(np.exp(2*mu + s*s)))]
    mean_analytic += [np.exp(mu+0.5*s*s)]
    sigma_analytic += [np.sqrt((np.exp(s**2)-1)*(np.exp(2*mu + s*s)))]

它给出了合理正确的平均值和标准差估计值，即使是在较高的西格玛值。在

关键是^{}使用MLE算法来估计参数。这完全不同于你最初的直接计算样本均值的方法。在

fit方法返回一个带有(sigma, loc=0, scale=exp(mu))的元组，这些元组是文档中指定的^{}对象的参数。在

我认为你应该调查一下你是如何估计均值和标准差的。分歧可能来自于算法的这一部分。在

其分歧可能有几个原因，至少考虑一下：

• 有偏估计量：你的估计量正确无偏吗？均值是无偏估计量（见下一节），但可能不是有效的
• 来自伪随机发生器的抽样异常值可能没有它们与理论分布进行比较时那么强烈：也许MLE比您的估计值更不敏感New MCVE bellow不支持这个假设，但是浮点运算错误可以解释为什么您的估计值被低估了
• 浮点运算错误新MCVE bellow强调这是您的问题的一部分。在

科学名言

似乎朴素的平均值估计（简单地取平均值），即使是无偏的，对于大西格玛来说，正确地估计平均值是低效的（见Qi Tang，Comparison of Different Methods for Estimating Log-normal Means，第11页）：

The naive estimator is easy to calculate and it is unbiased. However, this estimator can be inefficient when variance is large and sample size is small.

本文回顾了对数正态分布均值估计的几种方法，并以极大似然估计为例进行了比较。这就解释了为什么你的方法会随着西格玛的增加而产生漂移，而最大似然误差（MLE）就更好了——唉，对于大的N来说，这不是一个高效的方法。非常有趣的论文。在

统计因素

回顾：

• Lognormal是aheavy，长尾分布正偏态。一个结果是：随着形状参数sigma的增长，不对称性和skweness也会增加，异常值的强度也会增加。在
• 样本大小的影响：随着从一个分布中提取的样本数量的增长，期望有一个异常值的增加（范围也是如此）。在

建立一个新的MCVE

让我们建立一个新的MCVE使它更清晰。下面的代码从lognormal分布中绘制不同大小的样本（N范围在100和10000）中，其中形状参数变化（sigma范围介于0.1和10）之间，并且scale参数设置为酉的。在

import warnings
import numpy as np
from scipy import stats

# Make computation reproducible among batches:
np.random.seed(123456789)

# Parameters ranges:
sigmas = np.arange(0.1, 10.1, 0.1)
sizes = np.logspace(2, 5, 21, base=10).astype(int)

# Placeholders:
rv = np.empty((sigmas.size,), dtype=object)
xmean = np.full((3, sigmas.size, sizes.size), np.nan)
xstd = np.full((3, sigmas.size, sizes.size), np.nan)
xextent = np.full((2, sigmas.size, sizes.size), np.nan)
eps = np.finfo(np.float64).eps

# Iterate Shape Parameter:
for (i, s) in enumerate(sigmas):
    # Create Random Variable:
    rv[i] = stats.lognorm(s, loc=0, scale=1)
    # Iterate Sample Size:
    for (j, N) in enumerate(sizes):
        # Draw Samples:
        xs = rv[i].rvs(N)
        # Sample Extent:
        xextent[:,i,j] = [np.min(xs), np.max(xs)]
        # Check (max(x) + min(x)) <= max(x)
        if (xextent[0,i,j] + xextent[1,i,j]) - xextent[1,i,j] < eps:
            warnings.warn("Potential Float Arithmetic Errors: logN(mu=%.2f, sigma=%2f).sample(%d)" % (0, s, N))
        # Generate different Estimators:
        # Fit Parameters using MLE:
        fit = stats.lognorm.fit(xs, floc=0)
        xmean[0,i,j] = fit[2]
        xstd[0,i,j] = fit[0]
        # Naive (Bad Estimators because of Float Arithmetic Error):
        xmean[1,i,j] = np.mean(xs)*np.exp(-0.5*s**2)
        xstd[1,i,j] = np.sqrt(np.log(np.std(xs)**2*np.exp(-s**2)+1))
        # Log-transform:
        xmean[2,i,j] = np.exp(np.mean(np.log(xs)))
        xstd[2,i,j] = np.std(np.log(xs))

观察：新的MCVE在sigma > 4时开始发出警告。在

MLE作为参考

使用MLE估计形状和比例参数的效果很好：

以上两个数字显示：

• 估计误差随形状参数的增大而增大
• 估计误差随着样本量的增加而减小（CTL）

请注意，MLE也很适合形状参数：

浮点运算

这是值得画出的范围f绘制的样本与形状参数和样本大小：

或最小值和最大值之间的小数位数构成样本：

在我的设置中：

np.finfo(np.float64).precision  # 15
np.finfo(np.float64).eps        # 2.220446049250313e-16

这意味着我们最多可以处理15个有效数字，如果两个数字之间的幅度超过，那么最大的数字就会吸收较小的数字。在

一个基本的例子：如果我们只能保留四个有效数字，1 + 1e6的结果是什么？ 确切的结果是1,000,001.0，但必须四舍五入到1.000e6。这意味着：由于舍入精度，运算的结果等于最高的数字。它是有限精度算法的固有特性。在

上面的两个图结合统计考虑支持了您的观察，即增加N并不能改善对您的MCVE中sigma大值的估计。在

上面和下面的图显示了当sigma > 3我们没有足够的有效数字（小于5）来执行有效的计算。在

更进一步，我们可以说，估计值将被低估，因为最大数将吸收最小值，然后被低估的和将被N除以，使估计值在默认情况下有偏差。在

当形状参数足够大时，由于Arithmetic Float Errors，计算会有很大偏差。

它意味着使用以下数量：

np.mean(xs)
np.std(xs)

在计算时，由于xs中存储的值之间存在重要的差异，估计值将有很大的算术浮点误差。下图再现了您的问题：

如前所述，由于采样向量中的高值（少数异常值）会吸收较小的值（大多数采样值），所以估计值是默认的（不超过）。在

对数变换

如果我们应用logarithmic transformation，我们可以大大减少这种现象：

xmean[2,i,j] = np.exp(np.mean(np.log(xs)))
xstd[2,i,j] = np.std(np.log(xs))

然后对平均值is correct进行朴素估计，并且受算术浮点误差的影响要小得多，因为所有样本值都将在几十年内，而不是具有高于浮点运算精度的相对大小。在

实际上，对于每个N和sigma，使用log变换返回的平均值和std估计值的结果与MLE相同：

np.allclose(xmean[0,:,:], xmean[2,:,:])  # True
np.allclose(xstd[0,:,:], xstd[2,:,:])    # True

参考

如果你想在进行科学计算时寻找这类问题的完整和详细的解释，我建议你阅读read the excellent book：N.J.Higham，数值算法的准确性和稳定性，Siam，第二版，2002年。在

奖金

下面是图形生成代码的示例：

import matplotlib.pyplot as plt

fig, axe = plt.subplots()
idx = slice(None, None, 5)
axe.loglog(sigmas, xmean[0,:,idx])
axe.axhline(1, linestyle=':', color='k')
axe.set_title(r"MLE: $x \sim \log\mathcal{N}(\mu=0,\sigma)$")
axe.set_xlabel(r"Standard Deviation, $\sigma$")
axe.set_ylabel(r"Mean Estimation, $\hat{\mu}$")
axe.set_ylim([0.1,10])
lgd = axe.legend([r"$N = %d$" % s for s in sizes[idx]] + ['Exact'], bbox_to_anchor=(1,1), loc='upper left')
axe.grid(which='both')
fig.savefig('Lognorm_MLE_Emean_Sigma.png', dpi=120, bbox_extra_artists=(lgd,), bbox_inches='tight')