R 連續型資料機率分佈參數估計、適合度檢定教學與範例

使用 R 根據連續型的資料，估計機率分佈的參數，並且檢定模型的適合度，找出母體的機率分佈。

資料分析者在取得資料之後，通常都會想知道這些資料來自於什麼樣的機率分佈，在這篇教學中我們將使用 R 來解決這個問題。

產生測試資料

首先我們隨機產生一些 Gamma 分佈的資料，並加上一些常態分佈的雜訊，模擬測量上的隨機誤差。

# 設定亂數種子
set.seed(3)

# 產生的樣本數
N = 1000

# 產生 Gamma 分佈的樣本
x <- rgamma(N, shape = 10, scale = 3)

# 加上常態分佈的誤差
x <- x + rnorm(N, mean = 0, sd = 0.1)

處理步驟

一般來說，尋找資料機率分佈的步驟，大致上可以分為四個步驟：

視覺化分析，也就是畫出資料的直方圖，觀察資料的分佈狀況。
猜測該以何種機率分布來配適資料會比較好。
用各種統計的檢定方法，進行模型的適合度檢定。
若配適出來的模型不適合，則回到步驟 2，重新使用其他的模型配適看看。

視覺化分析

第一步的視覺化分析是最簡單的，只要將資料整理好，畫出直方圖即可：

# 載入 ggplot2 套件
library(ggplot2)

# 將資料整理成 Data Frame
my.df <- data.frame(x = x)

# 畫出直方圖
ggplot(my.df, aes(x)) + geom_histogram()

也可以畫出密度圖來觀看：

# 畫出密度圖
ggplot(my.df, aes(x)) + geom_density()

選擇機率分佈

機率分佈的選擇主要是靠著自己的經驗以及對於統計分佈的基礎知識來挑選，這裡由於我們的資料是自己產生的，所以我們知道資料的分佈是來自於 Gamma，形狀參數（shape parameter）為 10，尺度參數（scale parameter）為 3。

在實務上我們會先根據畫出來的圖形猜測資料的機率分佈，以這個例子來說資料有點像是常態分佈，但是資料的分佈沒有對稱，所以我們猜它可能是 Gamma 分佈。

決定分佈之後，再估計它的各個參數值，這時候可以嘗試使用 MASS 套件的 fitdistr 函數，使用最大概似（maximum-likelihood）的方式估計各個參數：

# 載入 MASS 套件
library(MASS)

# 以最大概似方法估計 Gamma 分佈的參數
fitdistr(x, dgamma, start = list(shape = 5, scale = 1))

     shape        scale   
  10.0295408    2.9545184 
 ( 0.4424277) ( 0.1336632)

估計出來的參數值也跟真正的值很接近。

適合度檢定

決定好機率分布之後，接著就要進行適合度檢定，檢驗看看這個猜測的機率分佈是否合理，而適合度檢定有非常多種，在 R 中也都有對應的套件可以使用，以下是各種適合度檢定的使用方式。

卡方檢定

卡方檢定是一個很常用的適合度檢定，不過它通常都是用於檢定離散型的分佈，若要拿來檢定連續型的資料分佈，必須先轉為離散型的分佈再進行卡方檢定：

# 將連續型的資料轉為離散型
brks <- seq(floor(min(x)), ceiling(max(x)))
x.brks.counts <- table(cut(x, breaks = brks, include.lowest=FALSE, right=FALSE))

# 產生機率分佈的理論值
library('zoo')
brks.cdf <- pgamma(brks, shape = 10, scale = 3)
null.probs <- rollapply(brks.cdf, 2, function(x) x[2]-x[1])

# 卡方檢定
chisq.test(x.brks.counts, p = null.probs, rescale.p = TRUE, simulate.p.value = TRUE)

        Chi-squared test for given probabilities with simulated p-value (based on
        2000 replicates)

data:  x.brks.counts
X-squared = 58.824, df = NA, p-value = 0.4428

Kolmogorov–Smirnov 檢定

Kolmogorov–Smirnov 檢定是一種簡單的無母數檢定方法，它會比較實際資料與指定分佈（或參考資料）的累積分布函數（CDF），判斷兩者是否有差異。

# Kolmogorov–Smirnov 檢定
ks.test(x, "pgamma", shape = 10, scale = 3)

        One-sample Kolmogorov-Smirnov test

data:  x
D = 0.023277, p-value = 0.6507
alternative hypothesis: two-sided

Cramér-von Mises 檢定

Cramér-von Mises 檢定也是一種用來檢定連續型機率分布的方法，在 R 中若要進行這種檢定，可以使用 goftest 套件。

# 載入 goftest 套件
library(goftest)

# Cramer-von Mises 檢定
cvm.test(x, "pgamma", shape = 10, scale = 3)

        Cramer-von Mises test of goodness-of-fit
        Null hypothesis: Gamma distribution
        with parameters shape = 10, scale = 3

data:  x
omega2 = 0.13188, p-value = 0.4499

Anderson-Darling 檢定

Anderson-Darling 檢定也是用來檢定資料是否來自於指定機率分佈的方法，在 R 中同樣可以使用 goftest 套件來處理。

# 載入 goftest 套件
library(goftest)

# Anderson-Darling 檢定
ad.test(x, "pgamma", shape = 10, scale = 3)

        Anderson-Darling test of goodness-of-fit
        Null hypothesis: Gamma distribution
        with parameters shape = 10, scale = 3

data:  x
An = 0.90978, p-value = 0.4082

關於 Anderson-Darling 與 Kolmogorov–Smirnov 兩種檢定的差異，可以參考 StackExchange 的文章。

測試各種機率分佈

通常猜測機率分佈與檢定的過程會重複好幾次，嘗試看看不同的機率分佈所配適的結果，以下是常態分佈、t 分佈與 Gamma 分佈的比較，表中的數值是檢定的 p-value。

	常態分佈	t 分佈	Gamma 分佈
參數估計	mean = 29.63 sd = 9.28	df = 4.66 ncp = 25.10	shape = 10.03 scale = 2.95
卡方檢定	0.001499	0.0004998	0.4518
Kolmogorov–Smirnov 檢定	0.05591	5.052e-08	0.3787
Cramér-von Mises 檢定	0.04791	6.341e-07	0.6286
Anderson-Darling 檢定	0.01805	6e-07	0.7933

比較過各種機率分佈與檢定的結果之後，就可以比較清楚看出來哪一個機率分佈是比較適合的了。

t 分佈四種檢定結果的 p-value 都非常低，所以非常不適合，而常態分佈的檢定情況則是剛好在邊緣地帶，只有一個 Kolmogorov–Smirnov 檢定勉強通過，而 Gamma 分佈的檢定結果則是非常好，所以 Gamma 還是最適合的機率分佈。

參考資料：Ensemble Blogging