統(tǒng)計學有兩大主要分支着绷，分別是描述性統(tǒng)計學和推斷統(tǒng)計學集乔。描述性統(tǒng)計學用于描述和概括數(shù)據(jù)的特征以及繪制各類統(tǒng)計圖表颇玷”颗總體數(shù)據(jù)，往往因為數(shù)據(jù)量太大而難以被獲取帖渠，所以就有了通過較小的樣本數(shù)據(jù)推測總體特性的推斷統(tǒng)計學谒亦。值得一提的是現(xiàn)今火熱的“大數(shù)據(jù)”一詞并不僅僅是指數(shù)據(jù)量大，在《大數(shù)據(jù)時代》一書中作者舍恩伯格強調(diào)“大數(shù)據(jù)”不是隨機樣本阿弃，而是所有數(shù)據(jù)诊霹，即總體羞延，這與傳統(tǒng)的統(tǒng)計研究方法是有很大區(qū)別的渣淳。

推斷統(tǒng)計學的一個研究方向就是用樣本數(shù)據(jù)估算總體的未知參數(shù)，稱之為參數(shù)估計伴箩。如果是用一個數(shù)值進行估計入愧，則稱為點估計；如果估計時給出的是一個很高可信度的區(qū)間范圍嗤谚，則稱為區(qū)間估計棺蛛。

本文先介紹了抽樣分布和中心極限定理，并用蒙特卡洛方法進行模擬巩步；然后引入置信區(qū)間的概念旁赊，并將之用于分析BRFSS數(shù)據(jù)中的BMI指數(shù)上。

首先依舊是導入相關Python模塊和數(shù)據(jù)椅野，其中brfss是專門用于讀取和清理美國行為風險因素監(jiān)控BRFSS調(diào)研數(shù)據(jù)的模塊终畅。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import brfss  # 該模塊用于處理BRFSS數(shù)據(jù)

%config InlineBackend.figure_format = 'retina'

df = brfss.ReadBrfss()  # 讀取BRFSS數(shù)據(jù)

這里主要關注反應胖瘦程度的BMI指數(shù)，并將這一數(shù)據(jù)存入bmi變量中竟闪，其數(shù)據(jù)量有40萬之多离福。

bmi = df.bmi.dropna()  # 取數(shù)據(jù)中的bmi列，并去除缺失值
len(bmi)  

405058

中心極限定理

如果我們將上述40萬多份的BMI數(shù)據(jù)看成是總體炼蛤，然后從中隨機抽取n個數(shù)據(jù)組成一份樣本妖爷，并計算該樣本的均值。重復這一過程1000次理朋，我們就得到了1000個樣本的均值分布絮识，即抽樣分布。

抽樣分布滿足中心極限定理嗽上，即在樣本量n越來越大時次舌，均值的抽樣分布將越來越接近正態(tài)分布，該分布的均值等于總體的均值炸裆；標準差垃它，在這里也稱為標準誤差SE滿足公式：

def sampling_distribution(data, sample_size=20, bins=40):
    '''抽樣分布模擬也殖，輸出均值、標準差以及直方圖务热、ECDF圖'''
    
    # 隨機抽樣
    sampling = [np.mean(np.random.choice(data, size=sample_size, replace=False)) for _ in range(1000)]  
    
    # 輸出總體和抽樣分布的均值忆嗜、標準差
    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):
    '''計算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，進行三次模擬冲甘。

sampling_distribution(bmi, sample_size=10)

mean of sample means: 27.95
population means: 28.04
Standard deviation of sample means: 2.04
Standard Error: 2.10

樣本量為10的抽樣分布

sampling_distribution(bmi, sample_size=20)

mean of sample means: 28.11
population means: 28.04
Standard deviation of sample means: 1.50
Standard Error: 1.49

樣本量為20的抽樣分布

sampling_distribution(bmi, sample_size=100)

mean of sample means: 28.05
population means: 28.04
Standard deviation of sample means: 0.69
Standard Error: 0.67

樣本量為100的抽樣分布

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

一般當n大于等于30時，樣本均值的抽樣分布近似為正態(tài)分布陶夜。此時我們可以用樣本的均值來估計總體的均值凛驮，這就是點估計的一種最簡單的方式。但從上述分布也可以看出条辟，樣本均值其實是以一定概率在總體均值附近浮動的黔夭，所以這就有了后面將要講的置信區(qū)間。

關于中心極限定理捂贿，還有一點需要強調(diào)的是纠修，無論變量原來的分布是什么樣的，其均值的抽樣分布在n足夠大時都會接近正態(tài)分布厂僧。比如我們研究BRFSS數(shù)據(jù)中人們每周運動的總時間（單位：分鐘）扣草，大部分人每周運動的時間少于500分鐘，而極少數(shù)人能達到3000分鐘颜屠，其直方圖反應數(shù)據(jù)大部分集中在左側(cè)辰妙，而右側(cè)有一條長長的尾巴。

exemin = df[df.exemin != 0].exemin.dropna()   # 提取鍛煉時間數(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時粗井，抽樣分布與正態(tà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)分布餐济，我們就來看看其均值和標準差的一些性質(zhì)。這里導入scipy的統(tǒng)計模塊胆剧，使用scipy.stats.norm()模擬標準正態(tài)分布絮姆，即均值為0，標準差為1秩霍。使用norm.pdf()計算概率密度篙悯，并繪制概率密度函數(shù)（PDF）圖。

import scipy.stats
norm = scipy.stats.norm()  # 標準正態(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()

標準正態(tài)分布

PDF圖中曲線下的面積代表了概率铃绒， 使用norm.cdf()可計算這部分面積鸽照，即累積概率分布。于是我們就可以得到變量距離均值在1個標準差范圍內(nèi)的概率為0.68匿垄，2個標準差范圍內(nèi)的概率是0.95移宅，3個標準差范圍內(nèi)的概率是0.997〈涣疲可見在正態(tài)分布中，數(shù)據(jù)主要集中在3個標準差之內(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

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

norm.interval(0.95)

(-1.959963984540054, 1.959963984540054)

norm.interval(0.99)

(-2.5758293035489004, 2.5758293035489004)

置信區(qū)間

在能夠計算正態(tài)分布中一定概率下對應的變量區(qū)間后席噩，我們再回到之前用樣本均值估計總體均值時遺留的問題班缰，即樣本的均值圍繞總體均值在一定范圍內(nèi)浮動的。我們需要估算總體均值在多大的概率下落在抽樣的隨機區(qū)間內(nèi)悼枢，這就是置信區(qū)間埠忘。

我們?nèi)匀粚?0多萬的bmi數(shù)據(jù)當成是總體，然后從中隨機抽取樣本量為100的數(shù)據(jù)馒索，根據(jù)中心極限定理繪制抽樣分布圖如下莹妒。

sample_size = 100    

# 計算總體的均值和標準差
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()

n=100抽樣分布

根據(jù)正態(tài)分布的性質(zhì)，在95%的概率下绰上，均值分布區(qū)間是(26.74, 29.35)旨怠。也就是說，在樣本量為100時蜈块，我們有95%的信心相信總體均值落在26.74和29.35之間鉴腻，這就是95%的置信區(qū)間迷扇。同理易核，99%的置信區(qū)間是(26.33, 29.76)骗村。注意這是在大樣本量的情況下，我們才能使用正態(tài)分布哑姚，而如果樣本量n小于30倦青，則需要采用t分布瓮床，此處就不展開了。

norm.interval(0.95)

(26.738141245959351, 29.346706751112283)

norm.interval(0.99)

(26.328305902131977, 29.756542094939658)

區(qū)間估計的應用

回到本系列文章一直在探索的一個問題产镐，即比較富人和普通人的BMI指數(shù)隘庄。此時，bmi數(shù)據(jù)不再當做總體看待癣亚，而是作為調(diào)查的樣本丑掺，總體是BRFSS數(shù)據(jù)針對的全體美國人。首先將bmi數(shù)據(jù)按照收入等級分為兩組述雾，即富人bmi數(shù)據(jù)和普通人bmi數(shù)據(jù)街州。

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

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

def mean_ci(data):
    '''給定樣本數(shù)據(jù)，計算均值95%的置信區(qū)間'''
    
    sample_size = len(data)
    std = np.std(data, ddof=1)  # 估算總體的標準差
    se = std / np.sqrt(sample_size)  # 計算標準誤差   
    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)匣掸。這兩個區(qū)間間隔的還比較遠趟紊，數(shù)值上差不多有1這么多。所以我們可以比較有信心的得出富人更瘦的結論碰酝。

mean_ci(bmi_rich)

(27.415906122294761, 27.485560606043915)

mean_ci(bmi_ord)

(28.509003170593907, 28.565637279855423)

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

mean_ci(bmi_rich[:500])

(27.849838839563304, 28.791561160436636)

mean_ci(bmi_ord[:500])

(28.200546441671069, 29.303493558328935)

此時富人bmi95%的置信區(qū)間為(27.85, 28.79)弄痹，而普通人bmi95%的置信區(qū)間為(28.20, 29.30)，很明顯這兩個區(qū)間都變大了嵌器。盡管富人的bmi指數(shù)仍有相對較小的趨勢肛真，但是這兩個區(qū)間有部分重合，這時我們就無法得出非乘剑肯定的結論了蚓让∏溃可見樣本量在做判斷時起著非常重要的作用，樣本越大历极，判斷越準確窄瘟，這也是與我們常識相符的。

小結

在這一篇中趟卸，我們了解了抽樣分布的概念蹄葱，中心極限定理的含義，正態(tài)分布的概率分布锄列，最重要的是使用置信區(qū)間的計算方法图云，通過樣本數(shù)據(jù)估算總體的均值范圍，至此我們進入了推斷統(tǒng)計學的領域邻邮。

針對富人是否更瘦這個問題上竣况，雖然使用了置信區(qū)間得出了較肯定的結論，但是仍然沒有對富人更瘦這個假設做出明確的判斷筒严。在下一篇中我們將會講到推斷統(tǒng)計學的另一個領域：假設檢驗丹泉，即對參數(shù)的假設值進行決策，屆時我們將和上述問題來個了斷鸭蛙。

本文代碼github

數(shù)據(jù)探索系列目錄：

1. 開篇：數(shù)據(jù)分析流程
2. 描述性統(tǒng)計分析
3. 統(tǒng)計分布
4. 參數(shù)估計（本文）
5. 假設檢驗

參考資料：

• 《Think Stats 2》
• 《統(tǒng)計學》摹恨，William Mendenhall著

致謝：
最后還是非常感謝解密大數(shù)據(jù)社群的小伙伴們給我的支持和鼓勵，讓我們一起成長规惰。

解密大數(shù)據(jù)社群公號