機率統計17——點估計和連續性修正

　　原文  |  https://mp.weixin.qq.com/s/NV3ThVwhM5dTIDQAWITSQQ

　　機率（probabilty）和統計（statistics）是兩個相近的概念，其實研究的問題恰好相反。

　　機率是使用一個已知參數的模型去預測這個模型所產生的結果，並研究結果的相關數字特徵，好比指望、方差等。假設如今已知一個射擊運動員的得分服從均值爲8.2，方差爲1.5的正態分佈，就能夠對這名運動員下一次射擊的得分狀況有個大體的預估。

　　統計與機率正好相反，是先有一堆數據，而後利用這些數據推測出模型和模型的參數。如今來了一名陌生的運動員，咱們對他一無所知，不過他宣稱本身是一名優秀的職業運動員。在進行了一系列射擊測試後，教練組收集到了一批這名運動員的成績，經過觀察數據，認爲他的成績符合正態分佈（也就是肯定了模型），以後再進一步經過數據推測模型參數的具體值。對於正態分佈來講，參數是均值和方差。

　　這樣看來，咱們碰到的大多數問題都是統計問題，雖然機率是隨機事件的客觀規律，但遺憾的是，這個規律老是做爲未知量出現的，做爲補償，咱們擁有一系列數據樣本，雖然這些樣本可能遠小於總體，但仍然能夠以點帶面，根據這些樣本對總體參數進行估計，得出關於總體機率分佈的近似值。這個根據樣本對整體參數進行估計的過程就是參數估計。根據參數性質不一樣，可分爲點估計和區間估計。

　　

　　點估計就是用樣本統計量的某一具體數值直接推斷未知的整體參數，獲得的是一個具體的參數值。咱們以前講過的矩估計、最大似然估計、貝葉斯估計都是點估計。

　　咱們以矩估計爲例，從新看看怎樣由樣本估計整體。

 

整體數量和樣本數量

　　咱們常常用n表示整體的數量，究竟什麼纔算整體呢？

　　整體老是給人以「多」的概念，但事實並不是如此，不一樣問題的「整體」數量可能相差很遠。好比一個啤酒廠一年生產的罐裝啤酒是1000萬，一個班級的學生數是60人，不管是1000萬仍是60，都是整體。用n表示整體的數量。

　　既然是用樣本估計整體，固然少不了抽樣，關於抽樣可參考：​數據分析(4)——閒話抽樣 | 看似公平的隨機抽樣是否真的公平？。 用m表示樣本的數量。

估計整體的均值

　　醫生的判斷很大程度依賴於血液檢測的結果。從抽血結束到取得報告單須要一段時間，這段時間並不肯定，也許運氣較好，十幾分鍾就拿到了，也可能等上一小時。如今獲得了一組樣本，X = {x₁, x₂, ……, x_m}，其中每一個數據表明了一名患者取得報告單前的等待時間。咱們的目標是根據這組樣本對總體的平均等待時間作一個估計。

　　計算方式很簡單，只須要計算樣本的均值就能夠了：

　　咱們認爲樣本的分佈與總體的分佈類似，這個均值是根據目前已知的數據對整體的最佳近似描述。這種用樣本均值估計整體均值的方法稱爲矩估計，估計的結果是整體均值的點估計量。

　　下面的代碼展示了總體分佈和抽樣分佈的關係：

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

fig = plt.figure(figsize=(10, 5))
plt.subplots_adjust(hspace=0.5)  # 調整子圖之間的上下邊距

mu, sigma_square = 30, 5 # 均值和方差
sigma = sigma_square ** 0.5 # 標準差
xs = np.arange(15, 45, 0.5)
ys = stats.norm.pdf(xs, mu, sigma)
ax = fig.add_subplot(2, 2, 1)
ax.plot(xs, ys, label='密度曲線')
ax.vlines(mu, 0, 0.2, linestyles='--', colors='r', label='均值')
ax.legend(loc='upper right')
ax.set_xlabel('X')
ax.set_ylabel('pdf')
ax.set_title('X~N($\mu$, $\sigma^2$), $\mu$={0}, $\sigma^2$={1}'.format(mu, sigma_square))

for i in [1, 2, 3]:
    m = 10 ** i # 樣本數量
    np.random.seed(m)
    X = stats.norm.rvs(loc=mu, scale=sigma, size=m) # 生成m個符合正態分佈的隨機變量
    X = np.trunc(X) # 對數據取整
    mu_x = X.mean() # 樣本均值
    ax = fig.add_subplot(2, 2, i + 1)
    ax.hist(X, bins=40)
    ax.set_xlabel('X')
    ax.set_ylabel('頻度')
    ax.set_title('m={0},樣本均值={1}'.format(m, mu_x))

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
# plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
plt.show()

　　能夠看到，抽取的樣本越多，樣本的分佈越接近總體的分佈。這裏有問題的是均值的符號，在各類資料中一會是μ，一會是戴帽子的μ，一會是x拔，到底用哪一個？

　　過去咱們一直說某些問題符合X~N(μ, σ²)的分佈，這裏μ是整體均值；在最大似然估計中得出的結果用表示，說明這是經過樣本對總體均值的一個點估計量。

估計整體的方差

　　假設咱們已經預先計算出了均值的點估計量，是否能夠像計算整體方差同樣計算樣本方差呢？

　　從上面的總體分佈圖能夠看到，大部分數據集中在均值附近，極端值出現的機率很低，這意味着對於抽樣來講，樣本數越小，抽到極端數值的可能性就越小。方差刻畫了數據相對於指望值的波動程度，因爲樣本的極端值出現的機率很低，所以樣本的波動極可能低於整體的波動，方差較總體方差更小。爲了應對這種狀況，咱們常常看到的是另外一個計算樣本方差的公式：

　　a/(m - 1)確定大於a/m，這使得②的結果稍大於①，m的值越小，①和②的差異越明顯。隨着樣本數量的增長，取得極端值的機會也變大，①和②的差異也會愈來愈小。將樣本的方差做爲整體方差的點估計量，一般用s²表示。

　　值得一提的是，若是咱們有m個樣本，在計算這些樣本的實際方差時，直接用①；若是是用這些樣本估計整體的方差，應該使用②。

估計整體的比例

　　不少人會在30分鐘以內取得報告單，一樣也有不少人要等更久。咱們能夠計算出樣本中成功人數（30分鐘以內拿到報告的人數）的比例，並用這個比例做爲整體機率的點估計量：

　　到目前爲止，點估計仍然很簡單，因此常常有人吐槽：機率這麼簡單的玩意有啥值得研究的？

樣本出現的機率

　　通過多年的統計分析，醫院已經明確告知，每一個患者都有50%的機率會在30分鐘內拿到報告單。咱們用p=50%表示整體中全部在30分鐘以內拿到報告的人數的佔比。若是把一個患者在30分鐘以內拿到報告看做成功，用隨機變量X表示m個樣本中的成功數量，那麼X符合參數爲m和p的二項分佈，X~B(m, p)，即成功次數符合試驗次數和成功率的二項分佈。保持試驗次數m不變，二項分佈近似於均值爲mp、方差爲mpq的正態分佈（q = 1 - p）。

　　下面的代碼畫出了二項分佈和其近似的正態分佈：

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

fig = plt.figure(figsize=(10, 6))
plt.subplots_adjust(hspace=0.8, wspace=0.3)  # 調整子圖之間的邊距

p = 0.5 # 每次試驗成功的機率
q = 1 - p # 每次試驗失敗的機率
m_list = [10, 15, 20] # 試驗次數
c_list = ['r', 'g', 'b'] # 曲線顏色
m_max = max(m_list)

# 二項分佈 X~B(m,p)
for i, m in enumerate(m_list):
    ax = fig.add_subplot(3, 2, i * 2 + 1)
    xs = np.arange(0, m + 1, 1) # 隨機變量的取值
    ys = stats.binom.pmf(xs, m, p) # 二項分佈 X~B(m,p)
    ax.vlines(xs, 0, ys, colors=c_list[i], label='m={}, p={}'.format(m, p))
    ax.set_xticks(list(range(0, m_max + 1, 2))) # 重置x軸座標
    ax.set_xlabel('X')
    ax.set_ylabel('pmf')
    ax.set_title('X~B(m, p)')
    ax.legend(loc='upper right')

# 保持二項分佈試驗的次數m不變，二項分佈近似於均值爲mp、方差爲mp(1-p)的正態分佈：
for i, m in enumerate(m_list):
    ax = fig.add_subplot(3, 2, i * 2 + 2)
    xs = np.arange(0, m + 1, 0.1) # 隨機變量的取值
    mu, sigma = m * p, (m * p * q) ** 0.5
    ys = stats.norm.pdf(xs, mu, sigma)
    ax.plot(xs, ys, c=c_list[i], label='m={}, p={}'.format(m, p))
    ax.set_xticks(list(range(0, m_max + 1, 2)))  # 重置x軸座標
    ax.set_xlabel('X')
    ax.set_ylabel('pdf')
    ax.set_title('X~N(mp, mpq)')
    ax.legend(loc='upper right')

plt.show()

　　某天來了20名患者，其中有12人在30分鐘以內拿到了報告單（12個成功）。根據二項分佈，這種狀況出現的機率是：

　　100天過去，天天都有20名患者接受驗血，x_i人在30分鐘內拿到了報告，天天的樣本都對應一個機率：

　　上式中全部m­_i的數量都是20，之因此用m­_i表示，是爲了強調雖然天天的樣本數量一致，但樣本自己是不一樣的。若是將這些機率也當作隨機變量，那麼這些變量也必然會符合某一個分佈，只要弄清這個分佈，就能回答產生某個樣本的機率。既然能夠經過樣本知道樣本中成功數量的佔比，那麼這個分佈也就等同於「樣本中成功數量的佔比」的機率。好比第10天的樣本中成功數量的佔比是p₁₀=45%，咱們的目標是瞭解p₁₀產生的機率有多大，即P(p₁₀)=？換句話說，咱們但願知道全部P_i(X=x_i)構成的分佈。

　　

　　咱們用p_s表示某個特定樣本中成功數量的佔比，藉助指望和方差來窺探p_s的分佈。一個明顯的關係是，若是整體中有50%的人能夠在30分鐘內拿到報告，那麼咱們也一樣指望在樣本中看到這個比例，這也是咱們可以用樣本估計整體的基礎。用隨機變量X表示樣本中成功的數量，p_s = X/m：

　　咱們已經知道X~B(m, p)，這裏m是樣本數量，p是每一個樣本成功的機率，是預先給出的。二項分佈的指望是E[x] = mp，方差是Var(X)=mpq，q = 1 – p，所以：

　　E[p_s]告訴咱們，樣本中成功數量的佔比與總體中成功數量的佔比一致；Var(p­_s)告訴咱們，m越大，p­_s的方差越小，樣本中成功數量的佔比越近整體中成功數量的佔比，用p_s來估計p越可靠。既然二項分佈X~B(m, p)能夠由X~N(mp, mpq)來近似，那麼p­_s =X/m也能夠由p­_s~N(p, pq/m)來近似。對於本例來講，p=0.25，pq/m=0.0125：

　　值得注意的是，比例的分佈刻畫的是樣本成功數佔比（即X/m）的變化狀況，而二項分佈刻畫的是特定數量的樣本中成功數（即X）的變化狀況。比例的取值範圍是[0, 1]，所以在描述p_s的分佈時，隨機變量的有效取值範圍是[0, 1]。當m固定時，每一個成功數佔比都表明一個特定的樣本，咱們能夠借用p_s的分佈計算從整體中抽樣出某個固定數量樣本的機率。

連續性修正

　　對於二項分佈來講，保持試驗次數n不變，二項分佈近似於均值爲np、方差爲npq的正態分佈。這裏特別強調了「近似於」，是由於二項分佈的隨機變量是離散型的，而正態分佈的隨機變量是連續型，可是這又有什麼關係呢？

　　這裏先要了解一下離散型分佈函數和連續型分佈函數的特色。對於連續型分佈來講，其分佈函數是用密度函數的積分表示的：

　　對於積分來講，a~b的區間與是否包含a點或b點沒什麼關係，對於連續型隨機變量的累積機率來講：

　　可是上式對於離散型隨機變量並不成立。下面是一個離散型分佈函數，縱座標的c.d.f是累積分佈函數（cumulative distribution function）的縮寫：

　　上圖向咱們展現了P(X < 1) = 0，P(X ≤ 1) = 0.5。這意味着對於離散型隨機變量來講，常常有P(x ≤ a) ≠ P(x < a)的狀況（並不老是不等，這要看a的取值，對於上圖來講，P(X < 1.5) = P(X ≤ 1.5)），而連續型隨機變量老是有P(x≤a) = P(x<a)。

　　

　　μ=50，σ²=25的正態分佈X~N(μ, σ²)能夠用來近似n=100，p=0.5的二項分佈X~B(n, p)，下圖是兩者的分佈函數（注意這裏的曲線是分佈函數，而不是密度函數）：

　　能夠看出，因爲二項分佈的離散型隨機變量只能取到整數，所以它的分佈函數是階梯狀的，而正態分佈的曲線穿過了每一個階梯的中心點，將階梯分紅了兩部分，左半部分離散分佈大於連續分佈，右半部分則相反：

　　分別用F_B(x)=P_B(X≤x)和F_n(x)=P_n(X≤x)表示二項分佈和正態分佈的分佈函數，對於整數x來講，在[x, x+0.5)區間內，F_B(x) > F_n(x)；在(x+0.5, x+1)區間內，F_B(x) < F_n(x)；只有在中心點，纔有F_B(x) = F_n(x)。

　　如今問題來了，用正態分佈去作近似的時候，若是直接用F_N(x)去近似P_B(X<x)，那麼結果會偏大；若是用F_N(x)去近似P_B(X≤x)，則結果會偏小：

　　時大時小並非個好主意，咱們想要的是一個一致的近似，要麼老是大，要麼老是小。一個辦法是對於X的正態連續性修正爲±0.5，即用F_N(x+0.5)去近似P_B(X < x)和P_B(X ≤ x)，獲得的結果不會偏小；或用F_N(x-0.5)去近似P_B(X < x)和P_B(X ≤ x)，獲得的結果不會偏大。這有點相似於用黎曼和計算積分時選用左矩形公式仍是右矩形公式：

　　回顧上一節的內容，咱們計算出了樣本佔比p­_s的正態分佈近似，p­_s的連續性修正爲：

　　

　　藉助連續性修正能夠求得：

　

---

　　出處：微信公衆號 "我是8位的"

　　本文以學習、研究和分享爲主，如需轉載，請聯繫本人，標明做者和出處，非商業用途！ 

　　掃描二維碼關注做者公衆號「我是8位的」