ICode9
精准搜索请尝试： 精确搜索
首页
编程语言>
数据库
系统相关
互联网
其他分享
Python
Java
JavaScript
android
PHP
首页 > 编程语言> 文章详细

python – 对数正态分布的均值和标准差与分析值不匹配

2019-05-27 06:42:29  阅读：709  来源： 互联网

标签：python floating-point numpy standard-deviation mean


作为我研究的一部分,我从对数正态分布中测量绘制的均值和标准差.给定基础正态分布的值,应该可以分析地预测这些量(如https://en.wikipedia.org/wiki/Log-normal_distribution所示).

但是,如下图所示,情况似乎并非如此.第一个图显示对数正态数据相对于高斯西格玛的平均值,而第二个图显示对数正态数据的西格玛与高斯的西格玛.显然,“计算”线非常显着地偏离“分析”线.

我将高斯分布的均值与μ= -0.5 * sigma ** 2的sigma相关联,因为这确保了对数正态场应该具有1的均值.注意,这是由物理学领域的推动我适用范围：例如,如果设置mu = 0.0,则仍会出现与分析值的偏差.

通过复制和粘贴问题底部的代码,应该可以重现下面的图表.任何关于可能导致这种情况的建议都将不胜感激！

高斯对数正态与西格玛的平均值：
Mean of lognormal vs sigma of gaussian

对数正态西格玛与高斯西格玛：
Sigma of lognormal vs sigma of gaussian

注意,为了产生上面的图,我使用了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

对数正态是正偏态和重尾分布.当对从高度偏斜的分布中抽取的样本执行浮点算术运算(例如sum,mean或std)时,采样向量包含的差异超过几个数量级(数十年).这使得计算不准确.

问题来自这两行：

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值,它也给出了合理正确的均值和标准差估计.

enter image description here
enter image description here

关键是fit使用MLE算法来估计参数.这与您直接执行样本均值的原始方法完全不同.

fit方法返回一个元组(sigma,loc = 0,scale = exp(mu)),它们是文档中指定的scipy.stats.lognorm对象的参数.

我认为您应该研究如何估算平均值和标准差.分歧可能来自算法的这一部分.

可能有几个原因导致它出现分歧,至少考虑：

>有偏估计：您的估算器是否正确且无偏见？均值是无偏估计(见下一节)但可能效率不高;
>伪随机生成器的采样异常值可能不像它们应该与理论分布相比那么强烈：可能MLE不如估计器敏感新MCVE波纹管不支持这个假设,但Float算术误差可以解释为什么你的估算器被低估了;
>浮动算术错误新的MCVE轰鸣声突出显示它是您问题的一部分.

科学的报价

似乎天真的平均估计量(简单地取均值),即使是无偏见的,也无法正确估计大型西格玛的均值(参见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更好地a对于大N来说它没有时间效率.非常有趣的论文.

统计考虑

回顾：

> Lognormal是heavy,长尾分布正向偏斜.一个结果是：随着形状参数sigma的增长,不对称性和skweness增长,异常值的强度也增加.
>样本大小的影响：随着从分布中抽取的样本数量增加,具有异常值的期望增加(程度也增加).

enter image description here
enter image description here

建立一个新的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开始发出警告. 4.

MLE作为参考

使用MLE估算形状和比例参数表现良好：

enter image description here
enter image description here

上面的两个数字显示：

>估计误差与形状参数一起增长;
>随着样本量的增加,估计误差减少(CTL);

注意比MLE还适合形状参数：

enter image description here

浮动算术

绘制样本的范围与形状参数和样本大小的关系是值得的：

enter image description here

或者,最小和最大数字之间的小数幅度构成样本：

enter image description here

在我的设置上：

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的估计.

上面和下面的数字表示比西格玛> 3我们没有足够的有效数字(小于5)来执行有效的计算.

更进一步,我们可以说估计量将被低估,因为最大数字将吸收最小值,然后低估的总和将除以N,使得估计量默认偏差.

当形状参数变得足够大时,由于Arithmetic Float Errors,计算被强烈偏置.

这意味着使用如下数量：

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

当计算估计将具有巨大的算术浮点错误,因为存储在xs中的值之间存在重要差异.下图重现了您的问题：

enter image description here
enter image description here

如上所述,估计是默认的(不过量),因为采样矢量中的高值(少数异常值)吸收小值(大多数采样值).

对数变换

如果我们应用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,采用对数变换返回与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')

标签：python,floating-point,numpy,standard-deviation,mean
来源： https://codeday.me/bug/20190527/1161587.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

关于我们 | 联系我们 | 留言反馈

专注分享技术，共同学习，共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有