作為我研究的一部分,我從對數(shù)正態(tài)分布中測量繪制的均值和標(biāo)準(zhǔn)差.給定基礎(chǔ)正態(tài)分布的值,應(yīng)該可以分析地預(yù)測這些量(如https://en./wiki/Log-normal_distribution所示).

但是,如下圖所示,情況似乎并非如此.第一個圖顯示對數(shù)正態(tài)數(shù)據(jù)相對于高斯西格瑪?shù)钠骄?而第二個圖顯示對數(shù)正態(tài)數(shù)據(jù)的西格瑪與高斯的西格瑪.顯然,“計算”線非常顯著地偏離“分析”線.

我將高斯分布的均值與μ= -0.5 * sigma ** 2的sigma相關(guān)聯(lián),因為這確保了對數(shù)正態(tài)場應(yīng)該具有1的均值.注意,這是由物理學(xué)領(lǐng)域的推動我適用范圍：例如,如果設(shè)置mu = 0.0,則仍會出現(xiàn)與分析值的偏差.

通過復(fù)制和粘貼問題底部的代碼,應(yīng)該可以重現(xiàn)下面的圖表.任何關(guān)于可能導(dǎo)致這種情況的建議都將不勝感激！

高斯對數(shù)正態(tài)與西格瑪?shù)钠骄担?br>

對數(shù)正態(tài)西格瑪與高斯西格瑪：

注意,為了產(chǎn)生上面的圖,我使用了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

對數(shù)正態(tài)是正偏態(tài)和重尾分布.當(dāng)對從高度偏斜的分布中抽取的樣本執(zhí)行浮點算術(shù)運算(例如sum,mean或std)時,采樣向量包含的差異超過幾個數(shù)量級(數(shù)十年).這使得計算不準(zhǔn)確.

問題來自這兩行：

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

因為ln包含非常低和非常高的值,其數(shù)量級遠(yuǎn)高于浮點精度.

可以使用以下謂詞輕松檢測到該問題以警告用戶其計算錯誤：

 (max(ln)   min(ln)) <= max(ln)

這在Strictly Positive Real中顯然是錯誤的,但在使用有限精度算術(shù)時必須考慮.

修改你的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值,它也給出了合理正確的均值和標(biāo)準(zhǔn)差估計.

關(guān)鍵是fit使用MLE算法來估計參數(shù).這與您直接執(zhí)行樣本均值的原始方法完全不同.

fit方法返回一個元組(sigma,loc = 0,scale = exp(mu)),它們是文檔中指定的scipy.stats.lognorm對象的參數(shù).

我認(rèn)為您應(yīng)該研究如何估算平均值和標(biāo)準(zhǔn)差.分歧可能來自算法的這一部分.

可能有幾個原因?qū)е滤霈F(xiàn)分歧,至少考慮：

>有偏估計：您的估算器是否正確且無偏見？均值是無偏估計(見下一節(jié))但可能效率不高;
>偽隨機生成器的采樣異常值可能不像它們應(yīng)該與理論分布相比那么強烈：可能MLE不如估計器敏感新MCVE波紋管不支持這個假設(shè),但Float算術(shù)誤差可以解釋為什么你的估算器被低估了;
>浮動算術(shù)錯誤新的MCVE轟鳴聲突出顯示它是您問題的一部分.

科學(xué)的報價

似乎天真的平均估計量(簡單地取均值),即使是無偏見的,也無法正確估計大型西格瑪?shù)木?參見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.

本文回顧了幾種估計對數(shù)正態(tài)分布均值的方法,并以MLE為參考進行比較.這就解釋了為什么你的方法隨著sigma的增加而變化,并且MLE更好地a對于大N來說它沒有時間效率.非常有趣的論文.

統(tǒng)計考慮

回顧：

> Lognormal是heavy,長尾分布正向偏斜.一個結(jié)果是：隨著形狀參數(shù)sigma的增長,不對稱性和skweness增長,異常值的強度也增加.
>樣本大小的影響：隨著從分布中抽取的樣本數(shù)量增加,具有異常值的期望增加(程度也增加).

建立一個新的MCVE

讓我們建立一個新的MCVE,使其更清晰.下面的代碼從對數(shù)正態(tài)分布中繪制不同大小的樣本(N范圍在100和10000之間),其中形狀參數(shù)變化(sigma范圍在0.1和10之間)并且比例參數(shù)被設(shè)置為單一的.

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=/).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))

觀察：當(dāng)西格瑪時,新的MCVE開始發(fā)出警告. 4.

MLE作為參考

使用MLE估算形狀和比例參數(shù)表現(xiàn)良好：

上面的兩個數(shù)字顯示：

>估計誤差與形狀參數(shù)一起增長;
>隨著樣本量的增加,估計誤差減少(CTL);

注意比MLE還適合形狀參數(shù)：

浮動算術(shù)

繪制樣本的范圍與形狀參數(shù)和樣本大小的關(guān)系是值得的：

或者,最小和最大數(shù)字之間的小數(shù)幅度構(gòu)成樣本：

在我的設(shè)置上：

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

這意味著我們最多可以使用15個有效數(shù)字,如果兩個數(shù)字之間的數(shù)量超過,則最大數(shù)字會吸收較小的數(shù)字.

一個基本的例子：如果我們只能保留四個有效數(shù)字,那么1 1e6的結(jié)果是什么？
確切的結(jié)果是1,000,001.0但必須四舍五入到1.000e6.這意味著：由于舍入精度,操作的結(jié)果等于最高數(shù).它是有限精度算術(shù)的固有特性.

上面的兩個數(shù)字結(jié)合統(tǒng)計學(xué)考慮支持您的觀察,即增加N不會改善MCVE中大值sigma的估計.

上面和下面的數(shù)字表示比西格瑪> 3我們沒有足夠的有效數(shù)字(小于5)來執(zhí)行有效的計算.

更進一步,我們可以說估計量將被低估,因為最大數(shù)字將吸收最小值,然后低估的總和將除以N,使得估計量默認(rèn)偏差.

當(dāng)形狀參數(shù)變得足夠大時,由于Arithmetic Float Errors,計算被強烈偏置.

這意味著使用如下數(shù)量：

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

當(dāng)計算估計將具有巨大的算術(shù)浮點錯誤,因為存儲在xs中的值之間存在重要差異.下圖重現(xiàn)了您的問題：

如上所述,估計是默認(rèn)的(不過量),因為采樣矢量中的高值(少數(shù)異常值)吸收小值(大多數(shù)采樣值).

對數(shù)變換

如果我們應(yīng)用logarithmic transformation,我們可以大大減少這種現(xiàn)象：

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

然后是平均值is correct的初始估計并且將更少受到算術(shù)浮點誤差的影響,因為所有樣本值都將在幾十年內(nèi),而不是具有高于浮點算術(shù)精度的相對幅度.

實際上,對于每個N和sigma,采用對數(shù)變換返回與MLE相同的均值和標(biāo)準(zhǔn)估計結(jié)果：

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

參考

如果您在進行科學(xué)計算時正在尋找關(guān)于此類問題的完整而詳細(xì)的解釋,我建議您參考read the excellent book：N.J. Higham,“數(shù)值算法的準(zhǔn)確性和穩(wěn)定性”,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')