• 如果您觉得本站非常有看点,那么赶紧使用Ctrl+D 收藏吧

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

开发技术 开发技术 7天前 4次浏览

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

目录
  • 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}^rfrac{(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)极大似然估计(hattheta),从而将分布变为具体的(f(x;hattheta)),再执行适当的分组让各个(p_i)已知,计算(chi^2)统计量。

这里的极大似然估计并不是样本的极大似然估计,而是对划分区间的极大似然估计,具体地,设(nu_j)为样本(X_1,cdots,X_n)中归类为第(j)类的个数,则似然函数为

[L=frac{n!}{nu_1!cdotsnu_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)}cdotfrac{partial p_i(theta)}{partialtheta}=0,
]

由此方程解出的(hattheta)是我们这里所需要的极大似然估计,它与由样本直接计算的极大似然估计是不一样的。

不过,需要注意的是,由于参数使用了估计量,所以(chi^2)统计量尽管计算方法相同,但却不再具有(chi^2_{r-1})的极限分布。Fisher指出:若存在参数空间(Theta)(s)维内点(theta),使得(X)的分布为(F(x;theta)),则对于(theta)的相合估计量(hattheta),在原假设成立的情况下,有

[K_nstackrel{d}to chi^2_{r-1-s}.
]

一般,矩估计不能用于待估参数的估计,寻找待估参数(hattheta)可以使用最小(chi^2)法,即定义

[K_n(theta)=sum_{i=1}^nfrac{(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_{1cdot} \
vdots & vdots & & vdots & & vdots & vdots \
i & n_{i1} & cdots & n_{ij} & cdots & n_{is} & n_{icdot} \
vdots & vdots & & vdots & & vdots & vdots \
r & n_{r1} & cdots & n_{rj} & cdots & n_{rs} & n_{rcdot} \
& n_{cdot 1} & cdots & n_{cdot j} & cdots & n_{cdot s} & n
end{matrix}
]

计算表明,应当用(n_{icdot}/n)来作为(p_i^{(1)})的估计值,(n_{cdot j}/n)来作为(p_j^{(2)})的估计值,于是计算所得的(chi^2)统计量为

[K_n=sum_{i=1}^rsum_{j=1}^sfrac{left(n_{ij}-nhat p_i^{(1)}hat p_j^{(2)}right)^2}{nhat p_i^{(1)}hat p_j^{(2)}}=nleft(sum_{i=1}^rsum_{j=1}^sfrac{n_{ij}^2}{n_{icdot}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_{1cdot}n_{2cdot}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_{icdot})事先选定的而非随机的。这一点关键的不同带来的问题是(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}^rsum_{j=1}^sfrac{(n_{ij}-n_{icdot}p_j^0)^2}{n_{icdot}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图检验,其原理是将要检验的样本的次序统计量,按照正态分布的分布函数反函数绘制到一张图上,如果这个散点图近似呈现一条直线,则认为这个样本服从正态分布。当然,这种方法有些主观,主要依靠人眼,不过由于十分直观,许多人也喜欢使用这个方法。

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


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

附录

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("ntPearson chi2 independence testnn")
  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("ntPearson chi2 consistency testnn")
  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")

程序员灯塔
转载请注明原文链接:数理统计15:拟合优度检验(2),列联表,正态性检验
喜欢 (0)