python - 对数正态分布的均值和标准差与分析值不匹配
问题描述
作为我研究的一部分,我从对数正态分布中测量平均数和标准差。给定基础正态分布的值,应该可以分析预测这些数量(如https://en.wikipedia.org/wiki/Log-normal_distribution中给出的)。
但是,如下图所示,情况似乎并非如此。第一个图显示了对数正态数据与高斯 sigma 的平均值,而第二个图显示了对数正态数据与高斯数据的 sigma。显然,“计算的”线与“分析的”线有很大的偏差。
我将高斯分布的平均值与 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;博士
对数正态分布是正偏态和重尾分布。在对从高度偏态分布中抽取的样本执行浮点算术运算(例如求和、均值或标准差)时,采样向量包含具有多个数量级(数十年)差异的值。这使得计算不准确。
问题来自这两行:
mean_calc += [np.average(ln)]
sigma_calc += [np.std(ln)]
因为ln
包含非常低和非常高的值,其数量级远高于浮点精度。
可以使用以下谓词轻松检测到该问题,以警告用户其计算错误:
(max(ln) + min(ln)) <= max(ln)
这在严格正实数中显然是错误的,但在使用有限精度算术时必须考虑。
修改你的 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
方法返回一个元组,(sigma, loc=0, scale=exp(mu))
它是scipy.stats.lognorm
文档中指定的对象的参数。
我认为您应该调查如何估计均值和标准差。分歧可能来自算法的这一部分。
它可能有几个不同的原因,至少考虑一下:
有偏估计:您的估计是否正确且无偏?均值是无偏估计量(见下一节),但可能效率不高;Pseudo Random Generator 的采样异常值可能不像理论分布那样强烈:也许 MLE 不如您的估计器敏感新的 MCVE 波纹管不支持这个假设,但 Float Arithmetic Error 可以解释为什么您的估计器被低估了;- 浮点算术错误下面的新 MCVE 突出显示它是您的问题的一部分。
科学报价
似乎朴素的均值估计器(简单地取均值),即使是无偏的,也无法正确估计大 sigma 的均值(参见Qi Tang,用于估计对数正态均值的不同方法的比较,第 11 页):
朴素估计量易于计算且无偏。但是,当方差较大且样本量较小时,此估计器可能效率低下。
论文回顾了几种估计对数正态分布均值的方法,并以MLE作为比较的参考。这解释了为什么您的方法会随着 sigma 的增加而漂移,而 MLE 会更好,因为它对于大 N 来说时间效率不高。非常有趣的论文。
统计考虑
回忆比:
- 对数正态分布是一个重且长尾的正偏态分布。一个后果是:随着形状参数 sigma 的增长,不对称性和偏度会增加,异常值的强度也会增加。
- 样本量的影响:随着从分布中抽取的样本数量的增加,对异常值的期望会增加(范围也会增加)。
构建新的 MCVE
让我们构建一个新的 MCVE 以使其更清晰。下面的代码从对数正态分布中抽取不同大小的样本(N
范围介于100
和之间),其中形状参数变化(范围介于和之间)并且比例参数设置为单一的。10000
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 也很适合 shape 参数:
浮点算术
值得绘制样本的范围与形状参数和样本大小的关系:
或者最小和最大数字之间的十进制大小形成样本:
在我的设置中:
np.finfo(np.float64).precision # 15
np.finfo(np.float64).eps # 2.220446049250313e-16
这意味着我们最多可以使用 15 个有效数字,如果两个数字之间的大小超过,那么最大的数字会吸收较小的数字。
一个基本的例子:1 + 1e6
如果我们只能保留四个有效数字,结果是什么?确切的结果是1,000,001.0
,但它必须四舍五入到1.000e6
。这意味着:由于舍入精度,运算结果等于最大数。它是有限精度算术固有的。
上面两个前两个数字结合统计考虑支持您的观察,即增加N
不会改善sigma
MCVE 中较大值的估计。
上图和下图显示了当sigma > 3
我们没有足够的有效数字(小于 5)来执行有效计算时。
此外,我们可以说估计器会被低估,因为最大的数字将吸收最小的数字,然后通过N
使估计器默认有偏差来划分被低估的总和。
当形状参数变得足够大时,由于算术浮点错误,计算会产生很大的偏差。
这意味着使用以下数量:
np.mean(xs)
np.std(xs)
由于存储在xs
. 下图重现了您的问题:
如前所述,估计是默认的(不过度),因为采样向量中的高值(少数异常值)吸收了小值(大部分采样值)。
对数变换
如果我们应用对数变换,我们可以大大减少这种现象:
xmean[2,i,j] = np.exp(np.mean(np.log(xs)))
xstd[2,i,j] = np.std(np.log(xs))
然后对平均值的天真估计是正确的,并且受算术浮点误差的影响要小得多,因为所有样本值都将位于几十年内,而不是具有高于浮点算术精度的相对幅度。
N
实际上,对于每个和,采用对数变换返回的均值和标准差估计结果与 MLE 相同sigma
:
np.allclose(xmean[0,:,:], xmean[2,:,:]) # True
np.allclose(xstd[0,:,:], xstd[2,:,:]) # True
参考
如果您在执行科学计算时正在寻找对此类问题的完整和详细的解释,我建议您阅读优秀的书籍:NJ Higham , Accuracy and Stability of Numerical Algorithms, Siam, Second Edition, 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')
推荐阅读
- javascript - 如何在反应中更新深层嵌套状态?
- netbeans - 在 NetBeans IDE 中使用 JDK11 启动 GlassFish 服务器
- python - 如何仅在 while 循环中运行 1 次操作
- javascript - 触发反应选择 onchange
- php - 在 php api 上验证图像文件时出错
- python - 安装 web3 时如何解决此错误
- java - 每次加载页面时,如何使用 Hibernate 将自动增量计数保存在数据库中?
- reactjs - @react-three/fiber 节点模块在 Expo 应用程序中不起作用
- css - 尝试设置输入类型范围的样式,但拇指未设置样式
- javascript - For 循环在 Firefox 上运行良好,但在 Chrome 上不行