數理統計15：拟合優度檢驗(2)，列聯表，正态性檢驗

數理統計15：拟合優度檢驗(2)，列聯表，正态性檢驗

本文我們繼續讨論拟合優度檢驗的相關問題。由于本系列為我獨自完成的，缺少審閱，如果有任何錯誤，歡迎在評論區中指出，謝謝！

目錄

• Part 1：分布未知的Pearson $\chi^2$檢驗
• Part 2：列聯表
  • Section 1：獨立性檢驗
  • Section 2：齊一性檢驗
• Part 3：正态性檢驗
• 附錄
  • 1、chisq.independence(mat)
  • 2、chisq.consistency(mat, prob)
  • 3、繪制QQ圖

Part 1：分布未知的Pearson \(\chi^2\)檢驗

在上篇文章中我們講到了分布已知的Pearson \(\chi^2\)檢驗，這是将觀測值分成\(r\)個組，根據實際頻數\(\nu_i\)和理論頻數\(np_i\)計算\(\chi^2\)統計量

\[K_n=\sum_{i=1}^r\frac{(\nu_i-np_i)^2}{np_i},

\]

将其近似為\(\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\)類的個數，則似然函數為

\[L=\frac{n!}{\nu_1!\cdots\nu_r!}(p_1(\theta))^{\nu_1}\cdots(p_r(\theta))^{\nu_r},

\]

對\(\ln L\)關于\(\theta\)求導（如果\(\theta\)是多元的，則是對向量求導），得到方程：

\[\sum_{i=1}^r \frac{\nu_i}{p_i(\theta)}\cdot\frac{\partial p_i(\theta)}{\partial\theta}=0,

\]

由此方程解出的\(\hat\theta\)是我們這裡所需要的極大似然估計，它與由樣本直接計算的極大似然估計是不一樣的。

不過，需要注意的是，由于參數使用了估計量，是以\(\chi^2\)統計量盡管計算方法相同，但卻不再具有\(\chi^2_{r-1}\)的極限分布。Fisher指出：若存在參數空間\(\Theta\)的\(s\)維内點\(\theta\)，使得\(X\)的分布為\(F(x;\theta)\)，則對于\(\theta\)的相合估計量\(\hat\theta\)，在原假設成立的情況下，有

\[K_n\stackrel{d}\to \chi^2_{r-1-s}.

\]

一般，矩估計不能用于待估參數的估計，尋找待估參數\(\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)}\)互相獨立。記

\[p_{ij}=\mathbb{P}(X^{(1)}=i,X^{(2)}=j),\quad \forall i=1,\cdots,r;j=1,\cdots,s.

\]

則\(X^{(1)},X^{(2)}\)獨立的充要條件是：存在\(p^{(1)}_1,\cdots,p_r^{(1)}\)和\(p_1^{(2)},\cdots,p_s^{(2)}\ge 0\)使得

\[\sum_{i=1}^r p_i^{(1)}=1,\quad \sum_{j=1}^s p_j^{(2)}=1,\\

p_{ij}=p_i^{(1)}\cdot p_j^{(2)},\quad \forall i=1,\cdots,r;j=1,\cdots,s.

\]

這樣便定義了一個含參數的拟合優度檢驗問題，完成了模型的轉化，接下來的推導步驟請盡可能了解，如果難以了解也可以直接記住結論。

對于列聯表而言，令\(n_{ij}\)為\(X^{(1)}=i,X^{(2)}=j\)的樣本個數，則可以作出這樣的列聯表：

\[\begin{matrix}

& 1 & \cdots & j & \cdots & s \\

1 & n_{11} & \cdots & n_{1j} & \cdots & n_{1s} & n_{1\cdot} \\

\vdots & \vdots & & \vdots & & \vdots & \vdots \\

i & n_{i1} & \cdots & n_{ij} & \cdots & n_{is} & n_{i\cdot} \\

\vdots & \vdots & & \vdots & & \vdots & \vdots \\

r & n_{r1} & \cdots & n_{rj} & \cdots & n_{rs} & n_{r\cdot} \\

& n_{\cdot 1} & \cdots & n_{\cdot j} & \cdots & n_{\cdot s} & n

\end{matrix}

\]

計算表明，應當用\(n_{i\cdot}/n\)來作為\(p_i^{(1)}\)的估計值，\(n_{\cdot j}/n\)來作為\(p_j^{(2)}\)的估計值，于是計算所得的\(\chi^2\)統計量為

\[K_n=\sum_{i=1}^r\sum_{j=1}^s\frac{\left(n_{ij}-n\hat p_i^{(1)}\hat p_j^{(2)}\right)^2}{n\hat p_i^{(1)}\hat p_j^{(2)}}=n\left(\sum_{i=1}^r\sum_{j=1}^s\frac{n_{ij}^2}{n_{i\cdot}n_{\cdot j}}-1 \right).

\]

在獨立性假設成立的情況下，\(K_n\)應當漸進服從***度為

\[rs-1-(r+s-2)=(r-1)(s-1)

\]

的\(\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\)的值：

\[K_n=\frac{n(n_{11}n_{22}-n_{12}n_{21})^2}{n_{1\cdot}n_{2\cdot}n_{\cdot 1}n_{\cdot 2}}.

\]

這是我們在高中時常用的計算式，而高中常用到的\(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\)個工廠産品品質相同，即

\[p_1(j)=\cdots =p_r(j),\quad \forall j=1,\cdots,s.

\]

齊一性檢驗與獨立性檢驗略有不同，其關鍵不同就在于此時\(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(j)=\cdots=p_r(j)=p_j^0,\quad \forall j=1,\cdots,s.

\]

此時，未知參數隻剩下\(p_{1}^{(1)},\cdots,p_{r}^{(1)}\)，因而極限分布的***度應該是\((rs-1)-(r-s)=r(s-1)\)，\(K_n\)的計算式也自然變成了

\[K_n=\sum_{i=1}^r\sum_{j=1}^s\frac{(n_{ij}-n_{i\cdot}p_j^0)^2}{n_{i\cdot}p_j^0}\stackrel{d}\to \chi^2_{r(s-1)}.

\]

類似編寫了

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")