統計學有兩大主要分支，分别是描述性統計學和推斷統計學。描述性統計學用于描述和概括數據的特征以及繪制各類統計圖表。總體數據，往往因為數據量太大而難以被獲取，所以就有了通過較小的樣本數據推測總體特性的推斷統計學。值得一提的是現今火熱的“大數據”一詞并不僅僅是指數據量大，在《大數據時代》一書中作者舍恩伯格強調“大數據”不是随機樣本，而是所有數據，即總體，這與傳統的統計研究方法是有很大區别的。

推斷統計學的一個研究方向就是用樣本數據估算總體的未知參數，稱之為參數估計。如果是用一個數值進行估計，則稱為點估計；如果估計時給出的是一個很高可信度的區間範圍，則稱為區間估計。

本文先介紹了抽樣分布和中心極限定理，并用蒙特卡洛方法進行模拟；然後引入置信區間的概念，并将之用于分析BRFSS數據中的bmi指數上。

首先依舊是導入相關Python模塊和數據，其中brfss是專門用于讀取和清理美國行為風險因素監控BRFSS調研數據的模塊。

import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import brfss # 該模塊用于處理BRFSS數據 %config InlineBackend.figure_format = 'retina' df = brfss.ReadBrfss() # 讀取BRFSS數據

這裡主要關注反應胖瘦程度的BMI指數，并将這一數據存入BMI變量中，其數據量有40萬之多。

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

405058

如果我們将上述40萬多份的BMI數據看成是總體，然後從中随機抽取n個數據組成一份樣本，并計算該樣本的均值。重複這一過程1000次，我們就得到了1000個樣本的均值分布，即抽樣分布。

抽樣分布滿足中心極限定理，即在樣本量n越來越大時，均值的抽樣分布将越來越接近正态分布，該分布的均值等于總體的均值；标準差，在這裡也稱為标準誤差SE滿足公式：

這裡使用蒙特卡洛模拟的方法，在40萬BMI數據中随機抽取n個數計算均值，并重複1000次，組成抽樣分布。以下的sampling_distribution()函數用于實現這一模拟過程，并繪制抽樣分布的直方圖和ECDF圖。

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的抽樣分布

觀察上面的輸出結果和圖形，我們發現随着樣本量的遞增，抽樣分布越來越靠近正态分布，其均值和标準差也越來越符合中心極限定理中給出的關系。

一般當n大于等于30時，樣本均值的抽樣分布近似為正态分布。此時我們可以用樣本的均值來估計總體的均值，這就是點估計的一種最簡單的方式。但從上述分布也可以看出，樣本均值其實是以一定概率在總體均值附近浮動的，所以這就有了後面将要講的置信區間。

關于中心極限定理，還有一點需要強調的是，無論變量原來的分布是什麼樣的，其均值的抽樣分布在n足夠大時都會接近正态分布。比如我們研究BRFSS數據中人們每周運動的總時間（單位：分鐘），大部分人每周運動的時間少于500分鐘，而極少數人能達到3000分鐘，其直方圖反應數據大部分集中在左側，而右側有一條長長的尾巴。

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

人們每周運動時間的分布

顯然這一數據分布并不滿足正态分布，但是我們采用上述相同的方法模拟其樣本均值的抽樣分布，在樣本量n為1000時，抽樣分布與正态分布符合的非常好。可見中心極限定理并不要求變量原來分布的樣子，這也正是其魅力所在。

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

運動時間均值的抽樣分布

既然中心極限定理中涉及了正态分布，我們就來看看其均值和标準差的一些性質。這裡導入scipy的統計模塊，使用scipy.stats.norm()模拟标準正态分布，即均值為0，标準差為1。使用norm.pdf()計算概率密度，并繪制概率密度函數（PDF）圖。

import scipy.stats norm = scipy.stats.norm() # 标準正态分布 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()可計算這部分面積，即累積概率分布。于是我們就可以得到變量距離均值在1個标準差範圍内的概率為0.68，2個标準差範圍内的概率是0.95，3個标準差範圍内的概率是0.997。可見在正态分布中，數據主要集中在3個标準差之内。

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

反過來，我們也可以通過概率來求變量分布的區間，這裡使用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)

在能夠計算正态分布中一定概率下對應的變量區間後，我們再回到之前用樣本均值估計總體均值時遺留的問題，即樣本的均值圍繞總體均值在一定範圍内浮動的。我們需要估算總體均值在多大的概率下落在抽樣的随機區間内，這就是置信區間。

我們仍然将40多萬的bmi數據當成是總體，然後從中随機抽取樣本量為100的數據，根據中心極限定理繪制抽樣分布圖如下。

sample_size = 100 # 計算總體的均值和标準差 mu = np.mean(bmi) se = np.std(bmi) / np.sqrt(sample_size) # 繪制正态分布的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抽樣分布

根據正态分布的性質，在95%的概率下，均值分布區間是(26.74, 29.35)。也就是說，在樣本量為100時，我們有95%的信心相信總體均值落在26.74和29.35之間，這就是95%的置信區間。同理，99%的置信區間是(26.33, 29.76)。注意這是在大樣本量的情況下，我們才能使用正态分布，而如果樣本量n小于30，則需要采用t分布，此處就不展開了。

norm.interval(0.95)

(26.738141245959351, 29.346706751112283)

norm.interval(0.99)

(26.328305902131977, 29.756542094939658)

回到本系列文章一直在探索的一個問題，即比較富人和普通人的BMI指數。此時，bmi數據不再當做總體看待，而是作為調查的樣本，總體是BRFSS數據針對的全體美國人。首先将bmi數據按照收入等級分為兩組，即富人bmi數據和普通人bmi數據。

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

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

def mean_ci(data): '''給定樣本數據，計算均值95%的置信區間''' 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%的置信區間為(27.42, 27.49), 普通人bmi95%的置信區間為(28.51, 28.57)。這兩個區間間隔得還比較遠，數值上差不多有1這麼多。所以我們可以比較有信心地得出富人更瘦的結論。

mean_ci(bmi_rich)

(27.415906122294761, 27.485560606043915)

mean_ci(bmi_ord)

(28.509003170593907, 28.565637279855423)

但要注意了，以上之所以能得到這麼肯定的結論，源于使用的樣本數據量非常大，這大大縮小了置信區間的範圍（這可以從中心極限定理中标準誤差的公式看出）。現在讓我們使用前500個數據，看看在樣本較少時會發生什麼情況。

mean_ci(bmi_rich[:500])

(27.849838839563304, 28.791561160436636)

mean_ci(bmi_ord[:500])

(28.200546441671069, 29.303493558328935)

此時富人bmi95%的置信區間為(27.85, 28.79)，而普通人bmi95%的置信區間為(28.20, 29.30)，很明顯這兩個區間都變大了。盡管富人的bmi指數仍有相對較小的趨勢，但是這兩個區間有部分重合，這時我們就無法得出非常肯定的結論了。可見樣本量在做判斷時起着非常重要的作用，樣本越大，判斷越準确，這也是與我們常識相符的。

在這一篇中，我們了解了抽樣分布的概念，中心極限定理的含義，正态分布的概率分布，最重要的是使用置信區間的計算方法，通過樣本數據估算總體的均值範圍，至此我們進入了推斷統計學的領域。

針對富人是否更瘦這個問題上，雖然使用了置信區間得出了較肯定的結論，但是仍然沒有對富人更瘦這個假設做出明确的判斷。在下一篇中我們将會講到推斷統計學的另一個領域：假設檢驗，即對參數的假設值進行決策，屆時我們将和上述問題來個了斷。

更多精彩资讯请关注tft每日頭條，我们将持续为您更新最新资讯!