在之前的內容中，我們完成了引數估計的步驟，今天起我們將進入假設檢驗部分，這部分內容可參照《數理統計學教程》（陳希孺、倪國熙）。由於本系列為我獨自完成的，缺少審閱，**如果有任何錯誤，歡迎在評論區中指出，謝謝**！ [TOC] ## Part 1：什麼是假設檢驗 假設檢驗是一種統計推斷方法，用來判斷樣本與樣本、樣本與總體的差異是由抽樣誤差引起還是本質差別造成的。其步驟，其實就是提出一個假設，然後用抽樣作為證據，判斷這個假設是正確的或是錯誤的，這裡判斷的依據就稱為該假設的一個檢驗。 假設檢驗在數理統計中有重要的用途，比如：橙子的平均重量是80斤，這就是一個假設。我們怎麼才能知道它是對的還是錯的？這需要我們對橙子總體進行抽樣，然後對樣本進行一定的處理，比如計算總體均值的區間估計，如果區間估計不包含80斤，就認為原假設不成立，便拒絕原假設。 當然，由於樣本具有隨機性，因此我們只是對該假設進行**檢驗**而不是證明，也就是說不論假設檢驗的結果是接受假設還是拒絕假設，都不能認為假設本身是正確的或是錯誤的。同時，假設的檢驗也不是唯一確定的，對任何假設都可以有無數種方案進行檢驗，比如上面的例子，$95\%$的區間估計是一種檢驗，$99\%$的區間估計也可以作為檢驗，$90\%$的當然也可以，只要事先確定了即可。 總之，要將實用問題轉化為統計假設檢驗問題處理，一般需要經歷以下幾個步驟： 1. 明確所要處理的問題，將其轉化為二元問題，只能用“是”和“否”來回答。 2. 設計適當的檢驗，規定假設的**拒絕域**，即拒絕假設時樣本$\boldsymbol{X}$會落入的區域範圍（當然也可以是統計量會落入的範圍，這兩個意思是一致的）。 3. 抽取樣本$\boldsymbol{X}$進行觀測，計算需要的統計量的值。 4. 根據樣本的具體值作出接受假設或者否定假設的決定。 以下是假設檢驗問題的一些常用概念： - **零假設**即**原假設**，指的是進行統計檢驗時預先建立的假設，一般是希望證明其錯誤的假設，用字母$H_0$表示。這種區分方式比較玄乎。 - **備擇假設**指的是與原假設所相對立的假設，如果我們拒絕了原假設，就意味著需要接受備擇假設。 - **拒絕域**指的是我們會在什麼時候拒絕原假設，也就是拒絕原假設時樣本所屬的集合，它是樣本空間的一個子集。 在實際生活中，可以作出的假設是多種多樣的，對不同的檢驗問題有不同的方法，我們接下來對其分別進行討論。 ## Part 2：擬合優度檢驗 **擬合優度檢驗**指的是如下的檢驗問題——對於總體$X$和已知的分佈$F$，$H_0$：$X$的分佈為$F$。即用$F(x)$這個分佈函式來擬合樣本$X_1,\cdots,X_n$，擬合的優良程度如何。 對於擬合優度檢驗，其檢驗思想是構造一個統計量$D$，這個統計量刻畫樣本$\boldsymbol{X}=(X_1,\cdots,X_n)$與理論分佈$F$的偏離程度，$D$越大代表理論分佈$F$與樣本表現出來的分佈越不相符。由於在$H_0$成立也就是$X$的分佈為$F$的條件下，$D$作為一個統計量具有確定的抽樣分佈，且這個統計量描述的是擬合的優良程度，自然$D$越小說明擬合得越好。因此，我們這樣定義**擬合優度**：在獲得$\boldsymbol{X}$的情況下計算出$D$的觀測值$D_0$，則擬合優度為 $$ p(D_0)=\mathbb{P}(D\ge D_0|H_0). $$ 直觀地想象擬合優度，由於作為偏離量$D$肯定是非負的，所以如果$D_0$越小，擬合優度就越大，同時這也代表著樣本與理論程度的偏離小，所以**擬合優度可以作為理論分佈$F$是$X$的分佈的可靠性度量**。有了擬合優度之後，如果擬合優度太小，我們就會否定$H_0$，否則接受$H_0$，這需要一個閾值$\alpha$。常常取$\alpha=0.05$，如果擬合優度小於$\alpha$就否定零假設，不過閾值的設定多少是有些主觀的，因此客觀的擬合優度往往比接受和否定的二元回答更有價值。 不過，我們現在還沒有討論偏離量具體應該如何定義，並且偏離量的定義方式也可以多種多樣，因此我們接下來將分情況討論擬合優度檢驗中，偏離統計量$D$的定義。 ### Section 1：離散型（分佈已知） 如果$X$是離散型隨機變數，則可以進一步區分為如下的兩種：有限與可列。我們先討論有限且分佈已知的情況。現在，原假設是$H_0$：$X$服從如下的分佈： $$ F:\begin{pmatrix} a_1 & a_2 & \cdots & a_r \\ p_1 & p_2 & \cdots & p_r \end{pmatrix},\quad \sum_{i=1}^r p_i=1. $$ 適用於這種情況的檢驗是**Pearson的$\chi^2$檢驗**，這是我們所接觸的第一個重要檢驗。以$\nu_i$記所有樣本中取值為$a_i$的個數，稱為**觀察頻數**，即 $$ \nu_i=\#\{j|X_j=a_i\}. $$ 如果樣本容量為$n$，則理論上，應當有以下關係式： $$ p_i\approx \frac{\nu_i}{n},\quad \nu_i\approx np_i. $$ 稱$np_i$為**理論頻數**，如果$H_0$成立，則$|np_i-\nu_i|$對所有$i$都不應該過大，但又不能簡單地取$(np_i-\nu_i)^2$，需要考慮權重問題。Pearson引進如下的統計量，常稱為**Pearson $\chi^2$統計量**，為 $$ K_n=\sum_{i=1}^r\frac{(np_i-\nu_i)^2}{np_i} $$ 這裡，使用$1/np_i$作為權重的調整值，構造出檢驗統計量。使用這個權重的原因是Pearson證明的定理：在$H_0$成立的條件下，有$K_n\stackrel{d}{\to }\chi^2_{r-1}$。由此，如果由樣本計算出$K_n$的觀測值是$k_0$，則擬合優度為 $$ p(k_0)=\mathbb{P}(K_n\ge k_0|H_0)\approx\mathbb{P}(\tilde\chi^2_{r-1}\ge k_0), $$ 這裡$\tilde\chi^2_{r-1}$指的是服從$\chi^2_{r-1}$分佈的隨機變數，這樣擬合優度就可以用$\chi^2$分佈的分位數表近似地計算。實際應用時，我們不一定需要這麼精確的擬合優度，故只要查詢閾值$\alpha$對應的分位數值$\chi^2_{\alpha}(r-1)$，並比較它與$k_0$的大小即可，當$\chi^2_{\alpha}(r-1)\le k_0$時拒絕$H_0$。 > 這裡，我們以書上的例題為例，介紹如何用R語言計算擬合優度。現在$n=556$， > $$ > \begin{matrix} > \nu_i \\ > p_i \\ > np_i > \end{matrix}\begin{pmatrix} > 315 & 108 & 101 & 32 \\ > 9/16 & 3/16 & 3/16 & 1/16 \\ > 312.75 & 104.25 & 104.25 & 34.75 > \end{pmatrix} > $$ > 呼叫自定義函式`pearson.chi2(obersved, prob)`，程式碼見附錄。這個函式將返回一個列表，其中我們所需要的是$\chi^2$統計量的值和擬合優度，故進行如下呼叫： > > ```R > > observed <- c(315, 108, 101, 32) > > prob <- c(9, 3, 3, 1) > > pearson.chi2(observed, prob) > > Pearson chi2 test > > The value of K: 0.470024 > p-value: 0.9254259 > Accept. > ``` > > 這裡，`The value of K`是卡方統計量的值，`p-value`是擬合優度。如果需要進一步使用返回值的其他資訊，這個函式會返回一個列表。 如果$X$是可列的離散隨機變數，則需要對觀測進行**分組**，將其劃分為$r$個區間，這樣就將其轉化為有限的情形。具體地，分組需要讓每一個組內的理論頻數和觀察頻數都不小於5為宜。 ### Section 2：連續型（分佈已知） 如果理論分佈是一個已知的連續型分佈，雖然分組法依然可以解決（具體可以參考書本內容），但是不論如何取分組間隔，都多少遺漏了資訊。對於連續型分佈的驗證，使用**柯爾莫哥洛夫檢驗（柯氏檢驗）**更為合適。為此，我們需要先介紹**經驗分佈函式**。 經驗分佈函式是一個結合了分佈函式的單調、非降、左連續特色的函式，不同的是它可以由樣本計算得出。簡而言之，經驗分佈函式就是將經驗分佈函式中，$F(x)=\mathbb{P}(X< x)$的概率變成頻率。 > 注意，蘇中根《概率論》中，分佈函式的定義是右連續的即$F(x)=\mathbb{P}(X\le x)$，而這裡經驗分佈函式的定義卻是左連續的。不過，將分佈函式定義成左連續或是右連續並不影響使用，我們不妨就遵照左連續的定義。 設$X_1,\cdots,X_n$是$X\sim F$中抽取的簡單隨機樣本，次序統計量為$(X_{(1)},\cdots,X_{(n)})$，則$\forall x$，定義以下函式為經驗分佈函式： $$ F_n(x)=\left\{\begin{array}l 0,& x\le X_{(1)};\\ \frac{k}{n},& X_{(k)}

對於每一個**確定的$x$**，易知 $$ F_n(x)=\frac{1}{n}\sum_{i=1}^n I_{-\infty0. $$ 注意這個$D_n$，它顯然可以用來描述經驗分佈與理論分佈的偏離程度，因此我們顯然可以用它入手來完善擬合優度檢驗理論。柯氏檢驗就是建立在經驗分佈函式基礎上的，此時我們取理論分佈為$F$計算$D_n=\sup\limits_{x\in\mathbb{R}}|F_n(x)-F(x)|$，再根據確定的閾值$\alpha$，讓$K(\lambda)=\alpha$找到合適的$\lambda$，則檢驗的臨界值就是 $$ D_{n,\alpha}=\frac{\lambda}{\sqrt{n}}=\frac{K^{-1}(1-\alpha)}{\sqrt{n}}. $$ 由於經驗分佈函式是階梯型的，真實的分佈函式又是單調的，所以$D_n$必定出在各個$X_i$的觀測值處，具體地有 $$ D_n=\max\left\{\left|F_0(x_{(i)})-\frac{i-1}{n} \right|,\left|F_0(x_{(i)})-\frac{i}{n} \right|,\quad \forall i=1,2,\cdots,n \right\}. $$ > R語言提供了Kolmogorov-Smirnov檢驗函式`ks.test(x, y)`，其中斯米爾諾夫檢驗適用於檢驗兩組樣本是否同分布的。以書上例題為例，可以如下呼叫： > > ```R > > num <- c(0.034, 0.437, 0.863, 0.964, 0.366, 0.469, 0.637, 0.632, 0.804, 0.261) # 書本誤將0.437打成了0.0437，可以在底下的計算表中看出 > > ks.test(num, qunif) > > One-sample Kolmogorov-Smirnov test > > data: num > D = 0.166, p-value = 0.9052 > alternative hypothesis: two-sided > ``` > > 如果我們需要驗證的已知分佈，則`ks.test`的第二個引數直接用分佈的名稱即可（採用預設引數，如果不是預設的$U(0,1)$或者$N(0,1)$則對樣本資料作變換）。 --- 今天，我們介紹了擬合優度檢驗的部分內容，由於篇幅所限，擬合優度檢驗的剩餘部分將放到下一篇文章中介紹。關於Pearson $\chi^2$獨立性檢驗的函式，我將其放到附錄中，在使用時可以用`source()`函式訪問，也可以直接將函式複製執行。 ## 附錄： ### 1、經驗分佈函式繪製 ```R rm(list=ls()) dev.off() opar <- par(no.readonly = T) x <- seq(0,6,0.001) y <- pexp(x, rate = 2) par(mfrow = c(2,2)) plot.ecdf(rexp(10, 2), cex=0.1, main = "10 samples") lines(x, y, col='red') plot.ecdf(rexp(50, 2), cex=0.1, main = "50 samples") lines(x, y, col='red') plot.ecdf(rexp(100, 2), cex=0.1, main = "100 samples") lines(x, y, col='red') plot.ecdf(rexp(200, 2), cex=0.1, main = "200 samples") lines(x, y, col='red') ``` ### 2、自定義函式`pearson.chi2()` ```R pearson.chi2 <- function(x, p, alpha=0.05){ # x是實際頻數，p是概率列表，alpha是閾值 rt <- list() r <- length(x) # 可能出現的情況 p <- p / sum(p) # 歸一化概率列表 n <- sum(x) # 總數 np <- n*p # 理論頻數 rt$observed <- x rt$dim <- r rt$total.observed <- n rt$prob <- p rt$expect <- np K <- 0 for (i in 1:r){ K = K + (np[i]-x[i])^2/np[i] } rt$pearson.chi2 <- K p.value <- 1 - pchisq(K, df=r-1) rt$p.value <- p.value cat("\n\tPearson chi2 test\n\n") cat("The value of K: ", K, "\n") cat("p-value: ", p.value, "\n") if (p.value > alpha){ cat("Accept.\n") } else{ cat("Reject.\n") } lst <- r