数理统计15:拟合优度检验(2),列联表,正态性检验

本文我们继续讨论拟合优度检验的相关问题。由于本系列为我独自完成的,缺少审阅,如果有任何错误,欢迎在评论区中指出,谢谢

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图检验,其原理是将要检验的样本的次序统计量,按照正态分布的分布函数反函数绘制到一张图上,如果这个散点图近似呈现一条直线,则认为这个样本服从正态分布。当然,这种方法有些主观,主要依靠人眼,不过由于十分直观,许多人也喜欢使用这个方法。

Rplot

拟合优度检验属于一种非参数假设检验,包括列联表中的相关问题。下一篇文章,我们将回到参数假设检验问题上,讨论正态分布参数的假设检验,这也是我们最常遇到的问题。

附录

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")
posted @ 2021-02-22 14:30  江景景景页  阅读(139)  评论(0编辑  收藏