數理統計15:擬合優度檢驗(2),列聯表,正態性檢驗
阿新 • • 發佈:2021-02-22
本文我們繼續討論擬合優度檢驗的相關問題。由於本系列為我獨自完成的,缺少審閱,**如果有任何錯誤,歡迎在評論區中指出,謝謝**!
[toc]
## 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()`函式進行獨立性檢驗,程式碼見附錄。
>
> ```R
> > 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$檢驗的例項。
```R
> 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()`函式:
```R
> 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)`
```R
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)`
```R
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圖
```R
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 = "re