本文我們繼續討論擬合優度檢驗的相關問題。由於本系列為我獨自完成的，缺少審閱，如果有任何錯誤，歡迎在評論區中指出，謝謝！
Part 1：分布未知的Pearson \(\chi^2\)檢驗
在上篇文章中我們講到了分布已知的Pearson \(\chi^2\)檢驗，這是將觀測值分成\(r\)個組，根據實際頻數\(\nu_i\)和理論頻數\(np_i\)計算\(\chi^2\)統計量
將其近似為\(\chi^2_{r-1}\)分布給出擬合優度。其條件是，每一個\(p_i\)都是已知的，如上例中孟德爾豌豆實驗的\(9:3:3:1\)，或者具體地服從參數為\(3\)的泊松分布，等等。
現在我們要討論的問題是，給定一定的樣本，要檢驗的問題是：它是否服從某一特定的分布族，如泊松分布族、均勻分布族、正態分布族，等等。它與上文中所討論的情況區別在於，此時不知道每一個具體的\(p_i\)是多少，自然就不能計算\(np_i\)是多少。
對於參數分布族而言，它們的未知性完全由少數幾個參數決定，這就給我們的擬合優度檢驗提供了思路：如果我們能把這些參數變成已知量，那么\(p_i\)就已知了，可以照常使用Pearson \(\chi^2\)檢驗。因此，如果給定了分布族但沒有給定具體的分布，則我們可以先估計出未知參數來。
這里使用的估計方法是點估計，我們選擇普適性更強的極大似然估計，對分布族\(f(x;\theta)\)，在給定樣本的情況下，可以先給出\(\theta\)的極大似然估計\(\hat\theta\)，從而將分布變為具體的\(f(x;\hat\theta)\)，再執行適當的分組讓各個\(p_i\)已知，計算\(\chi^2\)統計量。
這里的極大似然估計並不是樣本的極大似然估計，而是對划分區間的極大似然估計，具體地，設\(\nu_j\)為樣本\(X_1,\cdots,X_n\)中歸類為第\(j\)類的個數，則似然函數為
對\(\ln L\)關於\(\theta\)求導（如果\(\theta\)是多維的，則是對向量求導），得到方程：
由此方程解出的\(\hat\theta\)是我們這里所需要的極大似然估計，它與由樣本直接計算的極大似然估計是不一樣的。
不過，需要注意的是，由於參數使用了估計量，所以\(\chi^2\)統計量盡管計算方法相同，但卻不再具有\(\chi^2_{r-1}\)的極限分布。Fisher指出：若存在參數空間\(\Theta\)的\(s\)維內點\(\theta\)，使得\(X\)的分布為\(F(x;\theta)\)，則對於\(\theta\)的相合估計量\(\hat\theta\)，在原假設成立的情況下，有
一般，矩估計不能用於待估參數的估計，尋找待估參數\(\hat\theta\)可以使用最小\(\chi^2\)法，即定義
\[K_n(\theta)=\sum_{i=1}^n\frac{(\nu_i-np_i(\theta))^2}{np_i(\theta)}, \]尋找使得\(K_n(\theta)\)最小也就是偏離量最小的\(\theta\)作為參數估計，即求解
\[\frac{\partial K_n(\theta)}{\partial \theta}=0. \]這個方程組很難解。
使用極大似然法，將\(\theta\)關於划分區間的極大似然估計作為估計值會更簡單一些，並且也滿足條件。但是這個方程組也不是很容易求解。
為圖方便，我們可能會想着直接用樣本的極大似然估計（如正態分布中的\(\bar X\)和\(\frac{n-1}{n}S^2\)）來作為\(\theta\)的估計量，這樣計算量就小得多。但是此時\(K_n(\theta)\)不一定有極限分布\(\chi^2_{r-1-s}\)，更確切地說來，此時的擬合優度真值是介於\([\mathbb{P}(\chi^2_{r-s-1}\ge k_0),\mathbb{P}(\chi^2_{r-1}\ge k_0)]\)之間的。
關於這部分的內容，不再詳細展開，如果能夠計算出參數的估計值，就自然可以算出各個區間的概率，然后用pearson.chi2()函數即可計算。不過要注意，此時只有\(K_n\)的計算值是可用的，擬合優度需要用\(\chi^2_{r-s-1}\)的分布來計算。
Part 2：列聯表
列聯表指的是將每一個樣本的兩類指標\(A\)和\(B\)作統計，進而判斷兩個指標之間是否存在一定的關系，典型例子就是吸煙和肺癌之間的關聯。具體地，列聯表相關的問題又可以分為獨立性檢驗與齊一性檢驗，它們都可以使用\(\chi^2\)檢驗來完成，與擬合優度檢驗類似，又有所不同。
Section 1：獨立性檢驗
獨立性檢驗指的是，判斷某個個體的兩項指標\(A,B\)是否相關。一般地，可以將一個由大量個體構成的總體按照指標\(A\)划分成\(r\)級：\(A_1,\cdots,A_r\)，按照指標\(B\)划分成\(s\)級：\(B_1,\cdots,B_s\)，第\(i\)個個體的指標狀況為\((A_{r_i},B_{s_i})\)。
這將第\(i\)個個體看成隨機向量\(X_i\)，就有\(X_i=(r_i,s_i)\)。如果不同個體之間相互獨立，就可以將\(n\)個個體\(X_1,\cdots,X_n\)視為\(n\)個從總體\(X\)中抽取的簡單隨機樣本，指標\(A,B\)無關就意味着\(X=(X^{(1)},X^{(2)})\)的兩個分量\(X^{(1)},X^{(2)}\)相互獨立。記
則\(X^{(1)},X^{(2)}\)獨立的充要條件是：存在\(p^{(1)}_1,\cdots,p_r^{(1)}\)和\(p_1^{(2)},\cdots,p_s^{(2)}\ge 0\)使得
這樣便定義了一個含參數的擬合優度檢驗問題，完成了模型的轉化，接下來的推導步驟請盡可能理解，如果難以理解也可以直接記住結論。
對於列聯表而言，令\(n_{ij}\)為\(X^{(1)}=i,X^{(2)}=j\)的樣本個數，則可以作出這樣的列聯表：
計算表明，應當用\(n_{i\cdot}/n\)來作為\(p_i^{(1)}\)的估計值，\(n_{\cdot j}/n\)來作為\(p_j^{(2)}\)的估計值，於是計算所得的\(\chi^2\)統計量為
在獨立性假設成立的情況下，\(K_n\)應當漸進服從自由度為
的\(\chi^2\)分布。
R語言中有可以直接用於獨立性檢驗的函數，但是計算結果與我們實際的計算式略有不同，因此我們自編一個
chisq.independence()函數進行獨立性檢驗，代碼見附錄。> mat <- matrix(c(442, 38, 514, 6), nrow=2) > chisq.independence(mat) Pearson chi2 independence test The value of K: 27.13874 p-value: 1.893646e-07此時
p-value遠小於\(0.05\)，所以認為獨立性假設不成立。
特別在\(r=s=2\)時（考試可能出到的范圍），以下公式有助於簡化計算\(K_n\)的值：
這是我們在高中時常用的計算式，而高中常用到的\(3.841\)就是自由度為\(1\)，\(\alpha=0.05\)時的臨界值。
Section 2：齊一性檢驗
齊一性檢驗的提法是\(r\)個工廠\(A_1,\cdots,A_r\)生產的產品可以分為\(B_1,\cdots,B_s\)的\(s\)個等級，第\(i\)個工廠的\(j\)等品率為\(p_i(j)\)，要驗證的假設是\(r\)個工廠產品質量相同，即
齊一性檢驗與獨立性檢驗略有不同，其關鍵不同就在於此時\(A_1,\cdots,A_r\)的\(r\)個工廠中抽取的產品數不是隨機的，而是可以自由挑選的，也就是說\(n_{i\cdot}\)是事先選定的而非隨機的。這一點關鍵的不同帶來的問題是\(\chi^2_{(r-1)(s-1)}\)的極限分布是否依然存在，但幸好，有結論表明，極限分布依然是\(\chi^2_{(r-1)(s-1)}\)。對於這種情況，chisq.independence(mat)函數依然適用。
齊一性檢驗中還有一種情況，就是\(j\)等品的分布是已知的，即要驗證的假設增加為
此時，未知參數只剩下\(p_{1}^{(1)},\cdots,p_{r}^{(1)}\)，因而極限分布的自由度應該是\((rs-1)-(r-s)=r(s-1)\)，\(K_n\)的計算式也自然變成了
類似編寫了chisq.consistency(mat, prob)函數用於應對這種情況，不過由於書上沒有給出相關例題，我也難以保證這個函數的運行結果是\(100\%\)正確的。
Part 3：正態性檢驗
正態性檢驗是檢驗數據是否服從正態分布的，理論上Pearson \(\chi^2\)檢驗和柯氏檢驗都可以完成這個任務。我們將其單獨拿出來討論，是因為正態分布是一種重要的分布，故需要一些更有針對性的檢驗，而不是使用適用面廣的檢驗。
小樣本下，正態性檢驗的方法是\(W\)檢驗；大樣本下，正態性檢驗的方法是\(D\)檢驗。其原理我們不詳述，以下給出R語言中，如何使用這兩種檢驗方法進行正態性檢驗。
\(W\)檢驗又叫Shapiro-Wilk檢驗，在R中的函數為shapiro.test()。其原假設是\(H_0\)：\(X\)服從正態分布，計算出的\(W\)統計量越大，正態性越強，越應該接受原假設。盡管我們還沒有介紹什么是檢驗的p-value，不過讀者只需要知道，p-value越大越應該接受原假設，如果p-value小於給定的閾值\(\alpha\)（可以取\(0.05\)）就拒絕原假設。以下以書上的例子為例給出\(W\)檢驗的實例。
> v <- c(57, 66, 74, 77, 81, 87, 91, 95, 97, 109)
> shapiro.test(v)
Shapiro-Wilk normality test
data: v
W = 0.99134, p-value = 0.9982
\(D\)檢驗又叫D 'Agostino檢驗，在R語言的fBasics包中提供了dagoTest()函數：
> dagoTest(runif(100))
Test Results:
STATISTIC:
Chi2 | Omnibus: 26.2949
Z3 | Skewness: -0.3753
Z4 | Kurtosis: -5.1141
P VALUE:
Omnibus Test: 1.95e-06
Skewness Test: 0.7074
Kurtosis Test: 3.152e-07
可以看到，我們使用均勻分布隨機數時，綜合性檢驗、偏度檢驗、峰度檢驗有兩項都不能通過正態性檢驗。
還有一種直觀的檢驗方法：Q-Q圖檢驗，其原理是將要檢驗的樣本的次序統計量，按照正態分布的分布函數反函數繪制到一張圖上，如果這個散點圖近似呈現一條直線，則認為這個樣本服從正態分布。當然，這種方法有些主觀，主要依靠人眼，不過由於十分直觀，許多人也喜歡使用這個方法。
擬合優度檢驗屬於一種非參數假設檢驗，包括列聯表中的相關問題。下一篇文章，我們將回到參數假設檢驗問題上，討論正態分布參數的假設檢驗，這也是我們最常遇到的問題。
附錄
1、chisq.independence(mat)
chisq.independence <- function(mat){
rt <- list()
r <- nrow(mat)
s <- ncol(mat)
if (r == 1 || s == 1){
cat("錯誤：矩陣維度過低")
return(None)
}
rt$rows <- r
rt$cols <- s
sumrow <- apply(mat, 1, sum)
sumcol <- apply(mat, 2, sum)
n <- sum(mat)
rt$sumrow <- sumrow
rt$sumcol <- sumcol
rt$total <- n
K <- 0
for (i in 1:r){
for (j in 1:s){
K = K + mat[i, j]^2 / (sumrow[i] * sumcol[j])
}
}
K = n * (K - 1)
rt$X.square <- K
p.value <- 1 - pchisq(K, df=(r-1)*(s-1))
rt$p.value <- p.value
cat("\n\tPearson chi2 independence test\n\n")
cat("The value of K: ", K, "\n")
cat("p-value: ", p.value, "\n")
lst <- rt
}
2、chisq.consistency(mat, prob)
chisq.consistency <- function(mat, prob){
rt <- list()
r <- nrow(mat)
s <- ncol(mat)
if (s == 1){
cat("錯誤：矩陣維度過低")
return(None)
}
rt$rows <- r
rt$cols <- s
prob <- prob / sum(prob)
rt$prob <- prob
sumrow <- apply(mat, 1, sum)
sumcol <- apply(mat, 2, sum)
n <- sum(mat)
rt$sumrow <- sumrow
rt$sumcol <- sumcol
rt$total <- n
K <- 0
for (i in 1:r){
for (j in 1:s){
K = K + (mat[i, j] - sumrow[i] * prob[j])^2 / (sumrow[i] * prob[j])
}
}
rt$X.square <- K
p.value <- 1 - pchisq(K, df=r*(s-1))
rt$p.value <- p.value
cat("\n\tPearson chi2 consistency test\n\n")
cat("The value of K: ", K, "\n")
cat("p-value: ", p.value, "\n")
lst <- rt
}
3、繪制QQ圖
rm(list=ls())
dev.off()
opar <- par(no.readonly = T)
par(mfrow = c(2,2))
x1 <- rnorm(500, 3, 6)
qqnorm(x1, main = "Q-Q Graph of Norm")
qqline(x1, col = "red")
x2 <- runif(500, -3, 3)
qqnorm(x2, main = "Q-Q Graph of Unif")
qqline(x2, col = "red")
x3 <- rexp(500, 1)
qqnorm(x3, main = "Q-Q Graph of Exp")
qqline(x3, col = "red")
x4 <- rt(500, 5)
qqnorm(x4, main = "Q-Q Graph of T")
qqline(x4, col = "red")