拟合度检验富于教学内涵，没有应用意义。文中图解

拟合度检验中的正态随机向量，示以数值实例。文末给出（与全文其实无关的）应用建议和代码示例。

binom.test(x = 22,n = 32)
##
##  Exact binomial test
##
## data:  22 and 32
## number of successes = 22, number of trials = 32, p-value = 0.0501
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4999224 0.8388153
## sample estimates:
## probability of success
##                 0.6875

=32×0.5=16、标准误总体参数

=

。下图解释了为什么要做±0.5的校正——Pr(离散的二项随机数X≥22)≈Pr(对应连续正态随机数Z>21.5)，约等号左边的22恰好取到的概率近似于Pr(22.5≥Z>21.5)。

prop.test(x = 22,n = 32)
##
##  1-sample proportions test with continuity correction
##
## data:  22 out of 32, null probability 0.5
## X-squared = 3.7812, df = 1, p-value = 0.05183
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4986377 0.8325051
## sample estimates:
##      p
## 0.6875

=4.5000。

的右尾

prop.test(x = 22,n = 32,correct = F)
##
##  1-sample proportions test without continuity correction
##
## data:  22 out of 32, null probability 0.5
## X-squared = 4.5, df = 1, p-value = 0.03389
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5143332 0.8204746
## sample estimates:
##      p
## 0.6875

拟合度检验都不作连续性校正。

chisq.test(c(22,10))
##
##  Chi-squared test for given probabilities
##
## data:  c(22, 10)
## X-squared = 4.5, df = 1, p-value = 0.03389

拟合度检验

=0.5的

。

=1。缺省的虚无假设是三个百分比总体期望都是

。SPSS报告的

值与R的结果一致，都等于

。

chisq.test(c(17,10,5))
##
##  Chi-squared test for given probabilities
##
## data:  c(17, 10, 5)
## X-squared = 6.8125, df = 2, p-value = 0.03317

值应当与某一个二维正态向量的长度平方对应。为了使二水平情形与三水平情形更好类比，可以将二水平的问题表述为：Head的两种取值(0, 1)的百分比总体期望满足

=1。这正是下图斜率为-1的对角线上的点集。 一维随机向量以概率

出现蓝色箭头，以概率

出现粉色箭头。32枚个案有22次粉色箭头，10次蓝色箭头，平均值以黑点描出。黑点距离出发点

的标准误总体参数=无限枚箭头个案的总体标准差

。根据中心极限定理，计算出的从出发点到黑点向量以标准误作标准化，得到的向量长度平方值（在样本量不太小的时候）近似服从

分布。

取一般的数值比如0.4，以上的几何解读同样成立。

作图、验证的R代码如下。

mu1 <- 0.4; mu0 <- 1 - mu1;plot(c(0,1),c(1,0),col='grey',type='l',asp=1,xlab = "Head=1 Proportion",ylab="Head=0 Proportion")
abline(h=0,v=0,lty=3)
points(22/32,10/32,pch=19)points(mu1,mu0)
arrows(mu1,mu0,0,1,col='blue')
arrows(mu1,mu0,1,0,col='pink')chisq.test(xc <- c(22,10),p = mu<- c(mu1,mu0))
##
##  Chi-squared test for given probabilities
##
## data:  xc <- c(22, 10)
## X-squared = 11.021, df = 1, p-value = 0.0009009
sum((xc- sum(xc)*mu)^2 * (sum(xc)*mu)^-1)
## [1] 11.02083
(sigma2 <- mu1*sum((c(1,0)-mu)^2) + mu0*sum((c(0,1)-mu)^2))
## [1] 0.48
sum((xc/sum(xc) - mu)^2)/(sigma2/sum(xc))
## [1] 11.02083

,

,

)的取值空间。蓝、红、绿三种向量各自出现的概率为

。如果这三个总体参数都是1/3，三种向量的平均数根据中心极限定理近似服从（三维空间中）二维的正态分布。根据对称性，这个二维正态分布各向方差一致。可以计算：彩色个案的

=各彩色向量长度平方用概率加权求和，再除以自由度2。黑向量的标准误总体参数平方

=

。在样本量不太小的时候，黑向量长度平方/

得到的数值近似服从

分布。

## 数据与对称的参数
mu_neg <- 1/3; mu_neu <- 1/3;
mu_pos <- 1- mu_neg - mu_neu;
mu <- c(mu_neg,mu_neu,mu_pos)chisq.test(x = xc<-c(17,10,5),p = mu)
##
##  Chi-squared test for given probabilities
##
## data:  xc <- c(17, 10, 5)
## X-squared = 6.8125, df = 2, p-value = 0.03317
sum((xc- sum(xc)*mu)^2 * (sum(xc)*mu)^-1) # 6.8125
## [1] 6.8125## 三维作图
if(!require(plot3Drgl)){install.packages('plot3Drgl');require(plot3Drgl);}
## Loading required package: plot3Drgl
## Loading required package: rgl
## Loading required package: plot3D
open3d()
## wgl
##   1
plot3d(c(0,1,0,0),c(0,0,1,0),c(1,0,0,1),type='l',col='grey',asp=1,xlab='Tri<0 Prop',ylab = 'Tri=0 Prop',zlab = 'Tri>0 Prop')# points3d(mu_neg,mu_neu,mu_pos,pch=1)arrow3d(mu,c(1,0,0),col='blue')
arrow3d(mu,c(0,1,0),col='red')
arrow3d(mu,c(0,0,1),col='green')
arrow3d(mu,xc/sum(xc),col='black',type='r')## 对称情形的数值复核
(sigma2 <- sum(mu*c(sum((c(1,0,0)-mu)^2),sum((c(0,1,0)-mu)^2),sum((c(0,0,1)-mu)^2)))/2) # df=2
## [1] 0.3333333
sum((xc/sum(xc) - mu)^2)/(sigma2/sum(xc))
## [1] 6.8125

不对称，比如 (0.3, 0.5, 0.2)，三种向量的平均数根据中心极限定理仍然近似服从（三维空间中）二维的正态分布。但这个二维正态分布各向方差不再一致，需要在旋转后的新坐标架下估算图示中三角形截面上二维正态分布的方差-协方差矩阵总体参数。参考本专栏同主题文章

## 坐标旋转
(P_mt <- cbind(c(1,1,1)/ 3^.5,c(1,-1/2,-1/2)/ 1.5^.5,c(0,1,-1)/ 2^.5))
##           [,1]       [,2]       [,3]
## [1,] 0.5773503  0.8164966  0.0000000
## [2,] 0.5773503 -0.4082483  0.7071068
## [3,] 0.5773503 -0.4082483 -0.7071068
(prop_t <- ((xc <- c(17,10,5))/sum(xc)) %*% P_mt)
##           [,1]      [,2]      [,3]
## [1,] 0.5773503 0.2423974 0.1104854## 不对称的参数及其数值复核
mu_neg <- 0.3; mu_neu <- 0.5;
mu_pos <- 1- mu_neg - mu_neu;
mu <- c(mu_neg,mu_neu,mu_pos)chisq.test(x = xc,p = mu)
##
##  Chi-squared test for given probabilities
##
## data:  xc
## X-squared = 8.2604, df = 2, p-value = 0.01608
sum((xc- sum(xc)*mu)^2 * (sum(xc)*mu)^-1)
## [1] 8.260417
(mu_t <- mu %*% P_mt)
##           [,1]        [,2]     [,3]
## [1,] 0.5773503 -0.04082483 0.212132
(brg <- P_mt[,2:3]-cbind(rep(1,3))%*%mu_t[2:3])
##            [,1]       [,2]
## [1,]  0.8573214 -0.2121320
## [2,] -0.3674235  0.4949747
## [3,] -0.3674235 -0.9192388
(sigma2 <- t(brg)%*%diag(mu)%*%brg)
##             [,1]        [,2]
## [1,]  0.31500000 -0.07794229
## [2,] -0.07794229  0.30500000
rbind(prop_t[2:3] - mu_t[2:3])%*%solve(sigma2/sum(xc))%*%(prop_t[2:3] - mu_t[2:3])
##          [,1]
## [1,] 8.260417

t(sapply(list(c(1,2),c(1,3),c(2,3)),function(jy){Sigma2 <- t(brg2 <- (diag(3)-cbind(rep(1,3))%*%mu)[,jy]) %*%diag(mu) %*% brg2;z2 <- rbind(D2 <- (xc/sum(xc)-mu)[jy]) %*%solve(Sigma2/sum(xc)) %*% D2;c(jy, z2)})
)
##      [,1] [,2]     [,3]
## [1,]    1    2 8.260417
## [2,]    1    3 8.260417
## [3,]    2    3 8.260417

检验结果中，倒数第二行报告Exact Sig.为0.034，与多项分布检验缺省的报告结果(0.04119207≈0.0412)不一致。R中的EMT包提供了多项分布检验功能，与R原生的多项分布密度函数编程得到的结果一致。

chisq.test(xc<-c(17,10,5))
##
##  Chi-squared test for given probabilities
##
## data:  xc <- c(17, 10, 5)
## X-squared = 6.8125, df = 2, p-value = 0.03317if(!require(EMT)){  install.packages("EMT");require(EMT);}
## Loading required package: EMT
multinomial.test(xc ,prob =mu <-  rep(1/3,3))
##
##  Exact Multinomial Test, distance measure: p
##
##     Events    pObs    p.value
##        561   9e-04     0.0412## Verified (Exact Multinomial Test)
(d_cri <- dmultinom(xc,prob = mu))
## [1] 0.0009168089
dar <- array(dim=c(33,33))
for (i in 0:32) for (j in 0:32){k <- 32- i - j; dar[i+1,j+1]=ifelse(k<0,0,dmultinom(c(i,j,k),prob = mu))}
sum(dar) # all probabilities
## [1] 1
sum(dar[dar<=d_cri]) # p_value
## [1] 0.04119207

值定义，与多项分布密度定义的极端程度等级略有出入，正如《F 检验与Tukey HSD检验的备择假设方向图解》的图解，超出单维度之后，两种备择假设梯度不同。考虑到

值在学术史上被引入这个问题只是因为特定时代计算不便被迫查表，在多项分布检验各事件极端等级的比较中可能没必要越俎代庖替代多项分布密度本身。

值的概率，但只是多项分布的概率而不是

分布的概率。

## Exact Sig. of Chisq
multinomial.test(xc ,prob =mu,useChisq = T)
##
##  Exact Multinomial Test, distance measure: chisquare
##
##     Events    chi2Obs    p.value
##        561     6.8125     0.0345
(c_cri <- sum((xc-(ne <- mu*sum(xc)))^2 / ne))
## [1] 6.8125
dar.c <- darfor (i in 0:32) for (j in 0:32){k <- 32- i - j; dar.c[i+1,j+1]=ifelse(k<0,0,sum((c(i,j,k)-ne)^2 / ne));dar[i+1,j+1]=ifelse(k<0,0,dmultinom(c(i,j,k),prob = mu));}
sum(dar[dar.c>=c_cri]) # p_value
## [1] 0.03449689

alpha_familywise <- 0.05
xc<-c(17,10,5)
K <- length(xc)
alpha_corrected <- alpha_familywise/K
cl<-sapply(xc,FUN = function(x) binom.test(x,sum(xc),conf.level = 1-alpha_corrected)$conf)
rbind(Lower=cl[1,],Est=xc/sum(xc),Upper=cl[2,])
##            [,1]      [,2]       [,3]
## Lower 0.3128146 0.1367970 0.03996758
## Est   0.5312500 0.3125000 0.15625000
## Upper 0.7413646 0.5383827 0.36617910