參數(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ì)圖表。總體數(shù)據(jù)，往往因?yàn)閿?shù)據(jù)量太大而難以被獲取，所以就有了通過(guò)較小的樣本數(shù)據(jù)推測(cè)總體特性的推斷統(tǒng)計(jì)學(xué)。
推斷統(tǒng)計(jì)學(xué)的一個(gè)研究方向就是用樣本數(shù)據(jù)估算總體的未知參數(shù)，稱之為參數(shù)估計(jì)。如果是用一個(gè)數(shù)值進(jìn)行估計(jì)，則稱為點(diǎn)估計(jì)；如果估計(jì)時(shí)給出的是一個(gè)很高可信度的區(qū)間范圍，則稱為區(qū)間估計(jì)。
本文先介紹了抽樣分布和中心極限定理，并用蒙特卡洛方法進(jìn)行模擬；然后引入置信區(qū)間的概念，并將之用于分析BRFSS數(shù)據(jù)中的BMI指數(shù)上。
首先依舊是導(dǎo)入相關(guān)Python模塊和數(shù)據(jù)，順便看下數(shù)據(jù)量。
# 輸入
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import brfss
df = 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è)樣本的均值分布，即抽樣分布。
抽樣分布滿足中心極限定理，即在樣本量n越來(lái)越大時(shí)，均值的抽樣分布將越來(lái)越接近正態(tài)分布，該分布的均值等于總體的均值；標(biāo)準(zhǔn)差，在這里也稱為標(biāo)準(zhǔn)誤差SE滿足公式：
這里使用蒙特卡洛模擬的方法，在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.95
population means: 28.04
Standard deviation of sample means: 2.04
Standard Error: 2.10
mean of sample means: 28.11
population means: 28.04
Standard deviation of sample means: 1.50
Standard Error: 1.49
mean of sample means: 28.05
population means: 28.04
Standard deviation of sample means: 0.69
Standard 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ù)分布并不滿足正態(tài)分布，但是我們采用上述相同的方法模擬其樣本均值的抽樣分布，在樣本量n為1000時(shí)，抽樣分布與正態(tài)分布符合的非常好?？梢?jiàn)中心極限定理并不要求變量原來(lái)分布的樣子，這也正是其魅力所在。
# 輸入
sampling_distribution(exemin, sample_size=1000)
# 輸出
mean of sample means: 499.54
population means: 499.37
Standard deviation of sample means: 23.60
Standard 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.stats
norm = 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圖中曲線下的面積代表了概率， 使用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.683
2 sigma : 0.954
3 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)分布的PDF
norm = 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)誤差的公式看出）?，F(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í)相符的。
|
