專題

程式設計者的養成必須假以實戰，不能光寫或模仿一些小程式片段。本章設計幾個簡單但完整的專題，從問題出發，引導讀者從程式的角度來解決問題。需要用到的程式技巧已經在前面的章節介紹過，一方面複習學過的技巧與觀念，一方面逐漸熟悉解決問題的切入點與最終結果的表達（Presentation）。建議讀者看完題目後，先思考如何下手，甚至直接動手寫，以收實戰訓練之效。

專題 1： $\pi$ 的估計

單位圓的面積是 $\pi$ ，大家都知道 $\pi=3.14159...$ ，但是這個數字從何而來？一個簡單也不簡單的估計方法如下：

在單位圓周圍的外切正方形內，均勻的投擲 $N$ 個點，如下圖所示。
假設落在單位圓內的數量為 $M$ ，則依據面積比（單位圓 vs. 正方形面積）與落點分佈比的關係

$\frac{\pi}{4} \approx \frac{M}{N}$ 得出估計值 $\hat{\pi}=\frac{4M}{N}$ 。
顯然，當 $N$ 越大時（點愈多），估計值 $\hat{\pi}$ 將愈接近 $\pi$ 。請試著不斷調高 $N$ 值，觀察 $\hat{\pi}$ 的變化。
接著，選擇一個適當的、夠大的 $N$ ，做為估計 $\pi$ 的抽樣數，然後執行 $10000$ 次估計，結果表示為 $\hat{\pi}_i, i=1,2, \cdots, 10000$ 。計算這 $10000$ 個估計值的平均數與標準差，並繪製直方圖。

專題 2：雙樣本 T 檢定的 p-value 分佈

統計應用上經常需要檢定兩組資料的母體（來源）平均數是否相等，一般稱為雙樣本 T 檢定（Two-Sample T test），也是一般統計教科書的經典內容。這裡透過簡單的程式做實驗來理解其中的概念。以下是實驗的想法：

假設兩組資料都來自標準常態母體 $N(0,1)$ 且樣本數一樣同為 $n$ 。對這兩組資料進行雙樣本 T 檢定，將得到一個 p-value。想透過實驗知道這個 p-value 的分佈，換句話說，如果將這個 p-value 當作一個隨機變數，它的分配會是什麼？先猜猜看，再來寫程式做實驗。實驗設定 $n=10, 30, 100$ 。
承上，但假設兩組資料分別來自不同的常態分配 $N(0,1)$ 及 $N(0.5,1)$ （如下圖左）。在兩母體平均數不同的情況下，生成的資料在雙樣本 T 檢定下，p-value 將呈現什麼樣的分佈？
承上，繼續改善程式以便觀察當兩常態母體的平均數差距較大時（如下圖中、右），所生成的樣本在雙樣本 T 檢定下，p-value 分佈有什麼變化？
兩獨立樣本的 T 檢定可使用 scipy.stats.ttest_ind(x1, x2)
假設兩樣本分別來自常態 $N(0,1)$ 與 $N(0.5,1)$ ，且樣本數 $n=100$ ，計算指令 scipy.stats.ttest_ind 的檢定力。
同上，當樣本數 $n= 10, 20, 30, 50, 100, 300, 500$ 時，scipy.stats.ttest_ind 的檢定力分別是多少，繪製一張圖來展現樣本數與檢定力的關係。

專題 3：以蒙地卡羅實驗驗證 J-B 檢定統計量

令 $\{x_i, i=1,\cdots, n\}$ 代表來自標準常態 $N(0,1)$ 的 $n$ 個隨機樣本。統計量 $G_1$ 表示為

$G_1 = \sqrt{\frac{n}{6}} \hat{s}$ ，

其中 $\hat{s}$ 為偏態係數（skewness）的估計值（參考指令 scipy.stats.skew）。請利用蒙地卡羅模擬（Monte Carlo Simulation）驗證統計量 $G_1$ 服從標準常態 $N(0,1)$ 。其中蒙地卡羅模擬的環境設定（scenarios）為：
- 樣本數 $n = 10, 20, 30, 50, 100, 300, 500, 1000$ 。
- 針對每個樣本數 $n$ ，模擬次數皆為 $N=50,000$ 。
- 繪製 $n = 10$ 與 $n = 500$ 時，統計量 $G_1$ 的直方圖與 ECDF 圖。並分別畫上對應的標準常態 PDF 與 CDF 圖。
同上，但令統計量為

$G_2 = \sqrt{\frac{n}{24}} (\hat{k} - 3)$ ，

其中 $\hat{k}$ 為峰態係數（Kurtosis）的估計值（參考指令 scipy.stats.kurtosis）。同樣利用蒙地卡羅模擬，驗證統計量 $G_2$ 服從標準常態 $N(0,1)$ 。蒙地卡羅模擬的環境設定同上。
同上，但統計量為

$G_3 = G_1^2 + G_2^2 = \sqrt{\frac{n}{6}} \left(\hat{s}^2 +\frac{(\hat{k} - 3)^2}{4}\right)$ ，

同樣利用上述的蒙地卡羅模擬，驗證統計量 $G_3$ 服從卡方分配 $\chi^2(2)$ 。 $G_3$ 為著名的 J-B (Jarque-Bera) 常態檢定統計量。
將上述驗證程式改寫為一個副程式，假設取名為 stats, p_value = JB_test(x)，輸入參數 x 代表欲檢定是否為常態的一組資料。輸出兩個結果，stats 為 $G_3$ 檢定統計量的值，p_value 為檢定的 p-value。
接著檢驗檢定統計量 $G_3$ 的檢定力。採蒙地卡羅模擬方式，步驟如下：
- 從下列的分配母體中抽樣： $N(0,1), T(3), T(10), T(30), U(0,1), \chi^2(8)$ 。
- 抽樣數 $n=10, 20, 30, 50, 100, 300, 500$ 。
- 實驗次數 $N = 50000$ 。
- 型一誤 $\alpha = 0.05$ 。
- 對每個分配母體與樣本數，分別計算檢定力： $Power = P(Reject \; H_0 \;|\; H_a)$ ，其中 $H_0:$ 資料來自常態； $H_a:$ 資料來自其他分配。最後針對每個母體，繪製如下圖的 Power vs. sample size。觀察檢定力受樣本數與母體來源（與常態的相似度）的影響。其中 Y 軸必須選擇合適的範圍，方能呈現出清楚的 power 值。
- 其實，當 $H_a$ 來自常態時（也就是資料來自 $H_0$ 的意思），此時的 Power 又稱為顯著水準，且理論上， Power 應該維持在所設定的型一誤 $\alpha$ ，即 $0.05$ ，如下圖左。因此，當檢定統計量無法維持既定的顯著水準時，後續的檢定力也不用做了。因為檢定統計量是根據 $H_0$ 為真的條件下得到的，連「自己的出身」都維持不住，則檢定力也失去意義。
- 一般而言，檢定力會隨著樣本數增加而變大，亦即，樣本數大有利於辨認資料的「真偽」。

當資料來自 $H_0$ ，檢定力也稱「顯著水準」，理論上必須維持在 $\alpha$ 值

當資料來自 $H_0$ ，檢定力也稱「顯著水準」，理論上必須維持在 $\alpha$ 值

上述的 JB 統計量或稱 JB test 可以拿來與其他常態檢定的產品比試一下。譬如，scipy.stats 也提供了 Jarque–Bera test( jarque_bera), Kolmogorov-Smirnov test(kstest), Shapiro-Wilk test(shapiro) 與 Anderson-Darling test(anderson) 等常態檢定的指令。下圖以 Anderson-Darling test 與本專題的 JB test 比較（蒙地卡羅實驗 5 萬次），左圖為顯著水準的維持，右圖為對 T(3) 分配的檢定力。從左圖看出 scipy 的 anderson 指令在小樣本時，型一誤的維持較差，即便樣本數升高也比 JB test 差。另，右圖在對 T(3) 的檢定力也比 JB test 差，直到樣本數上升到 300 才趕上（本實驗沒有嘗試樣本數介於 100~300 之間）。雖然從下圖的比較中，JB test 優於 scipy 的 anderson，但這不意味著 JB test 比 Anderson-Darling test 好，很有可能是scipy 的這支 anderson 程式並沒有寫好。何況本專題的 JB test 在樣本數小於 2000 時，統計量並不滿足該有的卡方分配。讀者可以試試 scipy 的其他指令。
scipy 的這幾個指令都只能對一組資料進行檢定統計量與 p-val 的計算，並不能 broadcasting 到 5 萬次。因此執行效率較本專題的 JB test 差很多。讀者可以細心的觀察兩者間執行時間的差異。

汪群超 Chun-Chao Wang

專題

部落格統計

汪群超 Chun-Chao Wang

專題

Like this:

部落格統計