汪群超 Chun-Chao Wang

Dept. of Statistics, National Taipei University, Taiwan

專題

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

專題 1:\pi 的估計

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

  1. 在單位圓周圍的外切正方形內,均勻的投擲 N 個點,如下圖所示。

  2. 假設落在單位圓內的數量為 M,則依據面積比(單位圓 vs. 正方形面積)與落點分佈比的關係

    \frac{\pi}{4} \approx \frac{M}{N} 得出估計值 \hat{\pi}=\frac{4M}{N}

  3. 顯然,當 N 越大時(點愈多),估計值 \hat{\pi} 將愈接近 \pi。請試著不斷調高 N 值,觀察 \hat{\pi} 的變化。

  4. 接著,選擇一個適當的、夠大的 N,做為估計 \pi 的抽樣數,然後執行 10000 次估計,結果表示為 \hat{\pi}_i, i=1,2, \cdots, 10000。計算這 10000 個估計值的平均數與標準差,並繪製直方圖。


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

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

  1. 假設兩組資料都來自標準常態母體 N(0,1) 且樣本數一樣同為 n。 對這兩組資料進行雙樣本 T 檢定,將得到一個 p-value。想透過實驗知道這個 p-value 的分佈,換句話說,如果將這個 p-value 當作一個隨機變數,它的分配會是什麼?先猜猜看,再來寫程式做實驗。實驗設定 n=10, 30, 100

  2. 承上,但假設兩組資料分別來自不同的常態分配 N(0,1)N(0.5,1)(如下圖左)。在兩母體平均數不同的情況下,生成的資料在雙樣本 T 檢定下,p-value 將呈現什麼樣的分佈?

  3. 承上,繼續改善程式以便觀察當兩常態母體的平均數差距較大時(如下圖中、右),所生成的樣本在雙樣本 T 檢定下,p-value 分佈有什麼變化?

  4. 兩獨立樣本的 T 檢定可使用 scipy.stats.ttest_ind(x1, x2)

兩個不同母體間的差距示意圖

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

  1. \{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 = 10n = 500 時,統計量 G_1 的直方圖與 ECDF 圖。並分別畫上對應的標準常態 PDF 與 CDF 圖。

  2. 同上,但令統計量為

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

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

  3. 同上,但統計量為

    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) 常態檢定統計量。

  4. 將上述驗證程式改寫為一個副程式,假設取名為 stats, p_value = JB_test(x),輸入參數 x 代表欲檢定是否為常態的一組資料。 輸出兩個結果,stats 為 G_3 檢定統計量的值,p_value 為檢定的 p-value。

  5. 接著檢驗檢定統計量 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_a,檢定力隨著樣本數增加而變大
  • 上述的 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 差很多。讀者可以細心的觀察兩者間執行時間的差異。

%d bloggers like this: