---

參數(shù)估計(jì)

統(tǒng)計(jì)學(xué)有兩大主要分支，分別是描述性統(tǒng)計(jì)學(xué)和推斷統(tǒng)計(jì)學(xué)。描述性統(tǒng)計(jì)學(xué)用于描述和概括數(shù)據(jù)的特征以及繪制各類(lèi)統(tǒng)計(jì)圖表?？傮w數(shù)據(jù)，往往因?yàn)閿?shù)據(jù)量太大而難以被獲取，所以就有了通過(guò)較小的樣本數(shù)據(jù)推測(cè)總體特性的推斷統(tǒng)計(jì)學(xué)。

推斷統(tǒng)計(jì)學(xué)的一個(gè)研究方向就是用樣本數(shù)據(jù)估算總體的未知參數(shù)，稱(chēng)之為參數(shù)估計(jì)。如果是用一個(gè)數(shù)值進(jìn)行估計(jì)，則稱(chēng)為點(diǎn)估計(jì)；如果估計(jì)時(shí)給出的是一個(gè)很高可信度的區(qū)間范圍，則稱(chēng)為區(qū)間估計(jì)。

本文先介紹了抽樣分布和中心極限定理，并用蒙特卡洛方法進(jìn)行模擬；然后引入置信區(qū)間的概念，并將之用于分析BRFSS數(shù)據(jù)中的BMI指數(shù)上。
首先依舊是導(dǎo)入相關(guān)Python模塊和數(shù)據(jù)，順便看下數(shù)據(jù)量。

# 輸入import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as snsimport brfssdf = brfss.ReadBrfss()bmi = df.bmi.dropna() # 取數(shù)據(jù)中的bmi列，并去除缺失值print(len(bmi))# 輸出405058

中心極限定理

如果我們將上述40萬(wàn)多份的BMI數(shù)據(jù)看成是總體，然后從中隨機(jī)抽取n個(gè)數(shù)據(jù)組成一份樣本，并計(jì)算該樣本的均值。重復(fù)這一過(guò)程1000次，我們就得到了1000個(gè)樣本的均值分布，即抽樣分布。
抽樣分布滿(mǎn)足中心極限定理，即在樣本量n越來(lái)越大時(shí)，均值的抽樣分布將越來(lái)越接近正態(tài)分布，該分布的均值等于總體的均值；標(biāo)準(zhǔn)差，在這里也稱(chēng)為標(biāo)準(zhǔn)誤差SE滿(mǎn)足公式：

這里使用蒙特卡洛模擬的方法，在40萬(wàn)BMI數(shù)據(jù)中隨機(jī)抽取n個(gè)數(shù)計(jì)算均值，并重復(fù)1000次，組成抽樣分布。以下的sampling_distribution()函數(shù)用于實(shí)現(xiàn)這一模擬過(guò)程，并繪制抽樣分布的直方圖和ECDF圖。

def sampling_distribution(data, sample_size=20, bins=40):    #抽樣分布模擬，輸出均值、標(biāo)準(zhǔn)差以及直方圖、ECDF圖    # 隨機(jī)抽樣    sampling = [np.mean(np.random.choice(data, size=sample_size, replace=False)) for _ in range(1000)]      # 輸出總體和抽樣分布的均值、標(biāo)準(zhǔn)差    mu = np.mean(data)    se = np.std(data) / np.sqrt(sample_size)    print('mean of sample means: %.2f' % np.mean(sampling))    print('population means: %.2f' % mu)    print('Standard deviation of sample means: %.2f' % np.std(sampling))    print('Standard Error: %.2f' % se)    # 繪制抽樣分布的直方圖、ECDF圖    fig = plt.figure(figsize=(16,5))    p1 = fig.add_subplot(121)    plt.hist(sampling, bins=bins, rwidth=0.9)    plt.xlabel('sampling means')    plt.ylabel('counts')    p2 = fig.add_subplot(122)    plot_ecdf(sampling, xlabel='sampling means', label='sampling ')    sample = np.random.normal(mu, se, size=10000)    plot_ecdf(sample, xlabel='sampling means', label='normal distribution')    plt.show()def ecdf(data):    #計(jì)算ECDF    x = np.sort(data)    y = np.arange(1, len(x)+1) / len(x)    return (x,y)def plot_ecdf(data, xlabel=None , ylabel='ECDF', label=None):    #繪制ECDF圖    x, y = ecdf(data)    _ = plt.plot(x, y, marker='.', markersize=3, linestyle='none', label=label)    _ = plt.legend(markerscale=4)    _ = plt.xlabel(xlabel)    _ = plt.ylabel(ylabel)    plt.margins(0.02)

下面我們將樣本量n分別取為10、20、100，進(jìn)行三次模擬。

# 輸入sampling_distribution(bmi, sample_size=10)sampling_distribution(bmi, sample_size=20)sampling_distribution(bmi, sample_size=100)# 輸出mean of sample means: 27.95population means: 28.04Standard deviation of sample means: 2.04Standard Error: 2.10mean of sample means: 28.11population means: 28.04Standard deviation of sample means: 1.50Standard Error: 1.49mean of sample means: 28.05population means: 28.04Standard deviation of sample means: 0.69Standard Error: 0.67

觀察上面的輸出結(jié)果和圖形，我們發(fā)現(xiàn)隨著樣本量的遞增，抽樣分布越來(lái)越靠近正態(tài)分布，其均值和標(biāo)準(zhǔn)差也越來(lái)越符合中心極限定理中給出的關(guān)系。

一般當(dāng)n大于等于30時(shí)，樣本均值的抽樣分布近似為正態(tài)分布。此時(shí)我們可以用樣本的均值來(lái)估計(jì)總體的均值，這就是點(diǎn)估計(jì)的一種最簡(jiǎn)單的方式。但從上述分布也可以看出，樣本均值其實(shí)是以一定概率在總體均值附近浮動(dòng)的，所以這就有了后面將要講的置信區(qū)間。

關(guān)于中心極限定理，還有一點(diǎn)需要強(qiáng)調(diào)的是，無(wú)論變量原來(lái)的分布是什么樣的，其均值的抽樣分布在n足夠大時(shí)都會(huì)接近正態(tài)分布。比如我們研究BRFSS數(shù)據(jù)中人們每周運(yùn)動(dòng)的總時(shí)間（單位：分鐘），大部分人每周運(yùn)動(dòng)的時(shí)間少于500分鐘，而極少數(shù)人能達(dá)到3000分鐘，其直方圖反應(yīng)數(shù)據(jù)大部分集中在左側(cè)，而右側(cè)有一條長(zhǎng)長(zhǎng)的尾巴。

exemin = df[df.exemin != 0].exemin.dropna()   # 提取鍛煉時(shí)間數(shù)據(jù)，丟棄0或者缺失值plt.hist(exemin,bins=40, range=(0,3000), rwidth=0.9)  # 繪制直方圖plt.xlabel('exercise mins per week')plt.ylabel('counts')plt.show()

顯然這一數(shù)據(jù)分布并不滿(mǎn)足正態(tài)分布，但是我們采用上述相同的方法模擬其樣本均值的抽樣分布，在樣本量n為1000時(shí)，抽樣分布與正態(tài)分布符合的非常好。可見(jiàn)中心極限定理并不要求變量原來(lái)分布的樣子，這也正是其魅力所在。

# 輸入sampling_distribution(exemin, sample_size=1000)# 輸出mean of sample means: 499.54population means: 499.37Standard deviation of sample means: 23.60Standard Error: 23.75

正態(tài)分布的特性

既然中心極限定理中涉及了正態(tài)分布，我們就來(lái)看看其均值和標(biāo)準(zhǔn)差的一些性質(zhì)。這里導(dǎo)入scipy的統(tǒng)計(jì)模塊，使用scipy.stats.norm()模擬標(biāo)準(zhǔn)正態(tài)分布，即均值為0，標(biāo)準(zhǔn)差為1。使用norm.pdf()計(jì)算概率密度，并繪制概率密度函數(shù)（PDF）圖。

import scipy.statsnorm = scipy.stats.norm()  # 標(biāo)準(zhǔn)正態(tài)分布x = np.arange(-5, 5, 0.02)y = norm.pdf(x)  # 概率密度plt.plot(x,y)plt.axvline(x=0,ymax=0.95, linestyle='--', color='red', alpha=0.5)plt.axvline(x=1,ymax=0.59, linestyle='--', color='green')plt.axvline(x=-1,ymax=0.59, linestyle='--', color='green')plt.axvline(x=2,ymax=0.16, linestyle='--', color='blue')plt.axvline(x=-2,ymax=0.16, linestyle='--', color='blue')plt.margins(0.02)plt.show()

PDF圖中曲線(xiàn)下的面積代表了概率， 使用norm.cdf()可計(jì)算這部分面積，即累積概率分布。于是我們就可以得到變量距離均值在1個(gè)標(biāo)準(zhǔn)差范圍內(nèi)的概率為0.68，2個(gè)標(biāo)準(zhǔn)差范圍內(nèi)的概率是0.95，3個(gè)標(biāo)準(zhǔn)差范圍內(nèi)的概率是0.997?？梢?jiàn)在正態(tài)分布中，數(shù)據(jù)主要集中在3個(gè)標(biāo)準(zhǔn)差之內(nèi)。

# 輸入print('1 sigma : %.3f' % (norm.cdf(1) - norm.cdf(-1)))print('2 sigma : %.3f' % (norm.cdf(2) - norm.cdf(-2)))print('3 sigma : %.3f' % (norm.cdf(3) - norm.cdf(-3)))# 輸出1 sigma : 0.6832 sigma : 0.9543 sigma : 0.997

反過(guò)來(lái)，我們也可以通過(guò)概率來(lái)求變量分布的區(qū)間，這里使用norm.interval()，比如95%的情況下變量分布在-1.96到1.96之間，99%的情況下分布在-2.58到2.58之間。

# 輸入norm.interval(0.95)norm.interval(0.99)# 輸出(-1.959963984540054, 1.959963984540054)(-2.5758293035489004, 2.5758293035489004)

置信區(qū)間

在能夠計(jì)算正態(tài)分布中一定概率下對(duì)應(yīng)的變量區(qū)間后，我們?cè)倩氐街坝脴颖揪倒烙?jì)總體均值時(shí)遺留的問(wèn)題，即樣本的均值圍繞總體均值在一定范圍內(nèi)浮動(dòng)的。我們需要估算總體均值在多大的概率下落在抽樣的隨機(jī)區(qū)間內(nèi)，這就是置信區(qū)間。
我們?nèi)匀粚?0多萬(wàn)的bmi數(shù)據(jù)當(dāng)成是總體，然后從中隨機(jī)抽取樣本量為100的數(shù)據(jù)，根據(jù)中心極限定理繪制抽樣分布圖如下。

sample_size = 100    # 計(jì)算總體的均值和標(biāo)準(zhǔn)差mu = np.mean(bmi)se = np.std(bmi) / np.sqrt(sample_size)# 繪制正態(tài)分布的PDFnorm = scipy.stats.norm(mu, se)x = np.arange(26, 31, 0.01)y = norm.pdf(x)plt.plot(x,y)# 繪制抽樣分布的直方圖sample_size = 100    sampling = [np.mean(np.random.choice(bmi, size=sample_size, replace=False)) for _ in range(1000)]plt.hist(sampling, bins=40, rwidth=0.9, normed=True, alpha=0.7)plt.show()

根據(jù)正態(tài)分布的性質(zhì)，在95%的概率下，均值分布區(qū)間是(26.74, 29.35)。也就是說(shuō)，在樣本量為100時(shí)，我們有95%的信心相信總體均值落在26.74和29.35之間，這就是95%的置信區(qū)間。同理，99%的置信區(qū)間是(26.33, 29.76)。注意這是在大樣本量的情況下，我們才能使用正態(tài)分布，而如果樣本量n小于30，則需要采用t分布。

# 輸入norm.interval(0.95)norm.interval(0.99)# 輸出(26.738141245959351, 29.346706751112283)(26.328305902131977, 29.756542094939658)

區(qū)間估計(jì)的應(yīng)用

回到本系列文章一直在探索的一個(gè)問(wèn)題，即比較富人和普通人的BMI指數(shù)。此時(shí)，bmi數(shù)據(jù)不再當(dāng)做總體看待，而是作為調(diào)查的樣本，總體是BRFSS數(shù)據(jù)針對(duì)的全體美國(guó)人。首先將bmi數(shù)據(jù)按照收入等級(jí)分為兩組，即富人bmi數(shù)據(jù)和普通人bmi數(shù)據(jù)。

df2 = df[['bmi', 'income']].dropna()  # 提取數(shù)據(jù)中bmi和收入水平income這兩列，并忽略缺失值bmi_rich = df2[df2.income == 8].bmi   # 收入水平為8級(jí)的，認(rèn)為是富人bmi_ord = df2[df2.income != 8].bmi    # 收入水平為1-7級(jí)的，認(rèn)為是普通人群

以下定義了mean_ci()函數(shù)，根據(jù)置信區(qū)間的計(jì)算公式，計(jì)算95%置信度下均值所在的區(qū)間。

def mean_ci(data):    '''給定樣本數(shù)據(jù)，計(jì)算均值95%的置信區(qū)間'''    sample_size = len(data)    std = np.std(data, ddof=1)  # 估算總體的標(biāo)準(zhǔn)差    se = std / np.sqrt(sample_size)  # 計(jì)算標(biāo)準(zhǔn)誤差       point_estimate = np.mean(data)      z_score = scipy.stats.norm.isf(0.025)  # 置信度95%    confidence_interval = (point_estimate - z_score * se, point_estimate + z_score * se)    return confidence_interval

于是得到富人bmi95%的置信區(qū)間為(27.42, 27.49), 普通人bmi95%的置信區(qū)間為(28.51, 28.57)。這兩個(gè)區(qū)間間隔的還比較遠(yuǎn)，數(shù)值上差不多有1這么多。所以我們可以比較有信心的得出富人更瘦的結(jié)論。

# 輸入mean_ci(bmi_rich)mean_ci(bmi_ord)# 輸出(27.415906122294761, 27.485560606043915)(28.509003170593907, 28.565637279855423)

但要注意了，以上之所以能得到這么肯定的結(jié)論，源于使用的樣本數(shù)據(jù)量非常大，這大大縮小了置信區(qū)間的范圍（這可以從中心極限定理中標(biāo)準(zhǔn)誤差的公式看出）。現(xiàn)在讓我們使用前500個(gè)數(shù)據(jù)，看看在樣本較少時(shí)會(huì)發(fā)生什么情況。

# 輸入mean_ci(bmi_rich[:500])mean_ci(bmi_ord[:500])# 輸出(27.849838839563304, 28.791561160436636)(28.200546441671069, 29.303493558328935)

此時(shí)富人bmi95%的置信區(qū)間為(27.85, 28.79)，而普通人bmi95%的置信區(qū)間為(28.20, 29.30)，很明顯這兩個(gè)區(qū)間都變大了。盡管富人的bmi指數(shù)仍有相對(duì)較小的趨勢(shì)，但是這兩個(gè)區(qū)間有部分重合，這時(shí)我們就無(wú)法得出非?？隙ǖ慕Y(jié)論了。可見(jiàn)樣本量在做判斷時(shí)起著非常重要的作用，樣本越大，判斷越準(zhǔn)確，這也是與我們常識(shí)相符的。