对数正态分布的均值和标准差与分析值不匹配
Mean and standard deviation of lognormal distribution do not match analytic values
作为我研究的一部分，我测量了对数正态分布的均值和标准差。给定基础 正态 分布的值，应该可以分析预测这些量（如 https://en.wikipedia.org/wiki/Log-normal_distribution 处给出）。
然而，如下图所示，情况似乎并非如此。第一张图显示了对数正态数据的平均值与高斯西格玛的对比，而第二张图显示了对数正态数据与高斯西格玛的对比。显然，"calculated" 行与 "analytic" 行的偏差非常大。
我将高斯分布的均值与 sigma 相关联 mu = -0.5*sigma**2，因为这确保对数正态场的均值应为 1。请注意，这是受物理学领域的启发我在工作：如果您设置 mu=0.0 例如，与解析值的偏差仍然会发生。
通过复制并粘贴问题底部的代码，应该可以重现下面的情节。任何有关可能导致此问题的建议都将不胜感激！
对数正态平均值与高斯西格玛的对比：
对数正态西格玛与高斯西格玛：
注意，为了生成上面的图，我使用了 N=10000，但为了速度，在下面的代码中加入了 N=1000。
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()
TL;DR
对数正态分布是正偏态和重尾分布。当对从高度偏态分布中抽取的样本执行浮点算术运算（例如求和、均值或标准差）时，采样向量包含的值具有几个数量级（几十年）的差异。这使得计算不准确。
问题出在这两行：
mean_calc += [np.average(ln)]
sigma_calc += [np.std(ln)]
因为 ln 包含非常低和非常高的值，其数量级远高于浮点精度。
使用以下谓词可以很容易地检测到问题，以警告用户其计算错误：
(max(ln) + min(ln)) <= max(ln)
这在 Strictly Positive Real 中显然是错误的，但在使用有限精度算法时必须考虑。
修改您的 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)))]
即使对于高 sigma 值，它也能给出合理正确的均值和标准差估计。
关键是fit使用了MLE算法来估计参数。这与您直接执行样本平均值的原始方法完全不同。
fit 方法 returns 一个带有 (sigma, loc=0, scale=exp(mu)) 的元组，它们是文档中指定的 scipy.stats.lognorm 对象的参数。
我认为您应该研究如何估计均值和标准差。分歧可能来自你算法的这一部分。
造成分歧的原因可能有多种，至少考虑一下：
有偏估计器：你的估计器是否正确无偏？均值是无偏估计器（见下一节）但可能效率不高；与理论分布相比，伪随机生成器的样本异常值可能不如理论分布那么强烈：MLE 可能不如您的估计器敏感 新的 MCVE 波纹管不支持这个假设，但是 Float Arithmetic Error 可以解释为什么你的估计量被低估了；- Float 算术错误 New MCVE 波纹管强调这是您问题的一部分。
科学引述
似乎朴素的均值估计器（简单地取均值），即使是无偏的，也无法有效地正确估计大 sigma 的均值（参见 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作为参考进行比较。这解释了为什么你的方法随着 sigma 的增加而漂移，并且 MLE 坚持得更好，唉，对于大 N 来说，时间效率不高。非常有趣的论文。
统计考虑
回忆比：
- Lognormal is a heavy 和长尾分布呈正偏态。一个结果是：随着形状参数 sigma 的增加，不对称性和偏斜性增加，离群值的强度也增加。
- 样本大小的影响：随着从分布中抽取的样本数量的增加，出现异常值的期望值也会增加（范围也会增加）。
构建新的 MCVE
让我们构建一个新的 MCVE 以使其更清晰。下面的代码从形状参数变化的对数正态分布（sigma 范围在 0.1 和 10) 并且比例参数设置为单一的。
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 也很适合形状参数：
浮点运算
绘制样本的范围与形状参数和样本大小的关系图是值得的：
或样本中最小和最大数字之间的小数幅度：
关于我的设置：
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 并且受算术浮点误差的影响要小得多 因为所有样本值都在几十年内而不是相对幅度高于float 算术精度.
实际上，对于每个 N 和 sigma，采用对数变换 returns 均值和标准差估计的结果与 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')
作为我研究的一部分，我测量了对数正态分布的均值和标准差。给定基础 正态 分布的值，应该可以分析预测这些量（如 https://en.wikipedia.org/wiki/Log-normal_distribution 处给出）。
然而，如下图所示，情况似乎并非如此。第一张图显示了对数正态数据的平均值与高斯西格玛的对比，而第二张图显示了对数正态数据与高斯西格玛的对比。显然，"calculated" 行与 "analytic" 行的偏差非常大。
我将高斯分布的均值与 sigma 相关联 mu = -0.5*sigma**2，因为这确保对数正态场的均值应为 1。请注意，这是受物理学领域的启发我在工作：如果您设置 mu=0.0 例如，与解析值的偏差仍然会发生。
通过复制并粘贴问题底部的代码，应该可以重现下面的情节。任何有关可能导致此问题的建议都将不胜感激！
对数正态平均值与高斯西格玛的对比：
对数正态西格玛与高斯西格玛：
注意，为了生成上面的图，我使用了 N=10000，但为了速度，在下面的代码中加入了 N=1000。
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()
TL;DR
对数正态分布是正偏态和重尾分布。当对从高度偏态分布中抽取的样本执行浮点算术运算（例如求和、均值或标准差）时，采样向量包含的值具有几个数量级（几十年）的差异。这使得计算不准确。
问题出在这两行：
mean_calc += [np.average(ln)]
sigma_calc += [np.std(ln)]
因为 ln 包含非常低和非常高的值，其数量级远高于浮点精度。
使用以下谓词可以很容易地检测到问题，以警告用户其计算错误：
(max(ln) + min(ln)) <= max(ln)
这在 Strictly Positive Real 中显然是错误的，但在使用有限精度算法时必须考虑。
修改您的 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)))]
即使对于高 sigma 值，它也能给出合理正确的均值和标准差估计。
关键是fit使用了MLE算法来估计参数。这与您直接执行样本平均值的原始方法完全不同。
fit 方法 returns 一个带有 (sigma, loc=0, scale=exp(mu)) 的元组，它们是文档中指定的 scipy.stats.lognorm 对象的参数。
我认为您应该研究如何估计均值和标准差。分歧可能来自你算法的这一部分。
造成分歧的原因可能有多种，至少考虑一下：
有偏估计器：你的估计器是否正确无偏？均值是无偏估计器（见下一节）但可能效率不高；与理论分布相比，伪随机生成器的样本异常值可能不如理论分布那么强烈：MLE 可能不如您的估计器敏感新的 MCVE 波纹管不支持这个假设，但是 Float Arithmetic Error 可以解释为什么你的估计量被低估了；- Float 算术错误 New MCVE 波纹管强调这是您问题的一部分。
科学引述
似乎朴素的均值估计器（简单地取均值），即使是无偏的，也无法有效地正确估计大 sigma 的均值（参见 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作为参考进行比较。这解释了为什么你的方法随着 sigma 的增加而漂移，并且 MLE 坚持得更好，唉，对于大 N 来说，时间效率不高。非常有趣的论文。
统计考虑
回忆比：
- Lognormal is a heavy 和长尾分布呈正偏态。一个结果是：随着形状参数 sigma 的增加，不对称性和偏斜性增加，离群值的强度也增加。
- 样本大小的影响：随着从分布中抽取的样本数量的增加，出现异常值的期望值也会增加（范围也会增加）。
构建新的 MCVE
让我们构建一个新的 MCVE 以使其更清晰。下面的代码从形状参数变化的对数正态分布（sigma 范围在 0.1 和 10) 并且比例参数设置为单一的。
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 也很适合形状参数：
浮点运算
绘制样本的范围与形状参数和样本大小的关系图是值得的：
或样本中最小和最大数字之间的小数幅度：
关于我的设置：
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 并且受算术浮点误差的影响要小得多 因为所有样本值都在几十年内而不是相对幅度高于float 算术精度.
实际上，对于每个 N 和 sigma，采用对数变换 returns 均值和标准差估计的结果与 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')