本节书摘来自华章计算机《数学建模：基于r》一书中的第1章，第1.4节,作者：薛 毅 更多章节内容可以访问云栖社区“华章计算机”公众号查看。
本节要介绍的是另一类检验,其目标不是针对具体的参数,而是针对分布的类型.例如,通常假定总体分布具有正态性,而“总体分布为正态”这一断言本身在一定场合下就是可疑的,有待于检验.
假设根据某理论、学说甚至假定,某随机变量应当有分布f,现对x进行n次观察,得到一个样本x1,x2,…,xn,要据此检验h0:x具有理论分布f这里虽然没有明确指出对立假设,但可以说,对立假设为h1:x不具有理论分布f本问题的真实含义是估计实测数据与该理论或学说符合得怎么样,而不在于当认为不符合时,x可能备择的分布如何.故问题中不明确标出对立假设,反而使人感到提法更为贴近现实.
1.4.1 pearson拟合优度χ2检验
1.理论分布完全已知的情况
设x1,x2,…,xn为来自总体x的样本,将数轴(-∞,∞)分成m个区间i1=(-∞,a1), i2=［a1,a2), …, im=［am-1,∞)记这些区间的理论概率分别为p1, p2, …, pm, pi=p{x∈ii}, i=1,2,…,m记ni为x1,x2,…,xn中落在区间ii内的个数,则在原假设成立的条件下,ni的期望值为npi,ni与npi的差距(i=1,2,…,m)可视为理论与观察之间偏离的衡量,构造统计量k=∑mi=1(ni-npi)2npi（1.62）称k为pearson χ2统计量.pearson证明了,在原假设成立的条件下,当n→∞时,k依分布收敛于自由度为m-1的χ2分布.
在r中,使用chisq.test()函数计算pearson拟合优度χ2检验,其使用格式为chisq.test(x, y = null, correct = true,
p = rep(1/length(x), length(x)), rescale.p = false,
simulate.p.value = false, b = 2000)参数x为数值向量或矩阵(用于列联表检验),或者x和y同时为因子.y为数值向量,当x为矩阵时,y无效.如果x为因子,y必须为同样长度的因子.
correct为逻辑变量,表示是否对2×2列联表中的统计量作连续型修正,默认值为true.
p为向量,表示落在每个小区间内的理论概率,默认值为均匀分布.
rescale.p为逻辑变量,表示是否将输入的p标准化,默认值为false.当p的和不是1时,需要取rescale.p = true,此时重新计算p,使其和为1.
simulate.p.value为逻辑变量,表示是否使用monte carlo方法仿真计算p值,默认值为false.
b为正整数,表示monte carlo仿真的次数,默认值为2000.
例1.21 某消费者协会为了确定市场上消费者对5种品牌啤酒的喜好情况,随机抽取了1000名啤酒爱好者作为样本进行如下试验:5种啤酒按分别标注为a、b、c、d和e，让啤酒爱好者品尝，每人分别品尝过这5种啤酒后，选出一种作为他最喜欢的品牌.表1.8是根据样本资料整理得到的各种品牌啤酒爱好者的频数分布.试根据这些数据判断消费者对这5种品牌啤酒的爱好有无明显差异.
解 如果消费者对5种品牌啤酒喜好无显著差异,那么,就可以认为喜好这5种品牌啤酒的人呈均匀分布,即5种品牌啤酒爱好者人数各占20％.因此原假设为h0:喜好5种啤酒的人数服从均匀分布输入数据,调用chisq.test()函数计算.> x <- c(210, 312, 170, 85, 223)
chisq.test(x) chi-squared test for given probabilities
data: x
x-squared = 136.49, df = 4, p-value < 2.2e-16x-squared为χ2统计量,df为自由度,p-value为p值(=2.2×10-16)0.05,拒绝原假设,故认为消费者对5种品牌啤酒的喜好是有显著差异的.
例1.22 大麦的杂交后代关于芒性的比例应是无芒∶长芒∶短芒=9∶3∶4.实际观测值为335∶125∶160.试检验观测值是否符合理论假设.
解 根据题意h0:p1=916, p2=316, p3=416调用chisq.test()函数,其命令如下:> chisq.test(c(335, 125, 160), p = c(9, 3, 4)/16)
data: c(335, 125, 160)
x-squared = 1.362, df = 2, p-value = 0.5061p值(=0.5061)>0.05,接受原假设,即大麦的芒性比例符合9∶3∶4的比例.
例1.23 为研究电话总机在某段时间内接到的呼叫次数是否服从poisson分布,现收集了42个数据,如表1.9所示.通过对数据的分析,能否确认在某段时间内接到的呼叫次数服从poisson分布(α=0.1)？
解 编写相应的计算程序(程序名:exam0123.r)如下：x <- 0:6; y <- c(7, 10, 12, 8, 3, 2, 0)
f <- ppois(x, mean(rep(x, y))); m <- length(y)
p <- f［1］; p［m］ <- 1 - f［m-1］
for (i in 2:(m-1)) p［i］ <- f［i］ - f［i-1］
chisq.test(y, p=p)但计算结果会出现警告. chi-squared test for given probabilities
data: y
x-squared = 1.5057, df = 6, p-value = 0.9591
warning message:
chi-squared近似算法有可能不准in: chisq.test(y, p = p)为什么会出现这种情况呢?这是因为pearson拟合优度χ2检验要求:在分组后,每个组的理论频数和实际频数都要大于等于5,而后三组中出现的频数分别为3,2,0,均小于5.解决问题的方法是将后三组合成一组,此时的频数为5,满足要求.下面给出相应的r程序.z <- c(7, 10, 12, 8, 5) #% 重新分组
m <- length(z); p <- p［1:m-1］; p［m］ <- 1 - f［m-1］
chisq.test(z, p=p) #% 作检验计算得到p值(= 0.9696)0.1,因此,能确认在某段时间内接到的呼叫次数服从poisson分布. chi-squared test for given probabilities
data: z
x-squared = 0.5389, df = 4, p-value = 0.96962.理论分布依赖于若干个未知参数的情况
如果分布族f依赖于r个参数θ1,θ2,…,θr,要根据样本x1,x2,…,xn去检验假设h0:x的分布属于分布族{f(x;θ1,θ2,…,θr)}解决这个问题的步骤是,先通过样本作出(θ1,θ2,…,θr)的极大似然估计(θ1,θ2,…,θr)再检验假设h0:x具有分布f(x;θ1,θ2,…,θr)然后再按理论分布已知的情况进行处理,所不同的是由式(1.62)得到的统计量k服从自由度为m-1-r的χ2分布,即自由度减少了r.
从这种角度来讲,例1.23的计算是有问题的,因为在计算过程中,是用样本均值代替总体的数学期望,因此,需要减少一个自由度.
为了便于此类问题的计算,这里编写一个小程序(程序名:chi2gof.r)来完成这项工作.chi2gof <- function(chi2test, nparams){
df <- chi2testparameter - nparams;
pval <- 1 - pchisq(chi2teststatistic, df)
names(pval) <- "p-val"
c(chi2teststatistic, pval, df)
}在程序中, chi2test为chisq.test()函数生成的对象,nparams为未知参数的个数.函数的返回值有三个,分别是χ2统计量、p值和自由度.
用chi2gof()函数再计算例1.23,程序和计算结果如下：> source("chi2gof.r")
chi2gof(chisq.test(z, p = p), nparams = 1)
x-squared p-val df
0.5388945 0.9102670 3.0000000此时,自由度为3而不是4, 虽然p值有所减少但结论不变.
1.4.2 kolmogorov-smirnov检验
kolmogorov(科尔莫戈罗夫)拟合优度检验是kolmogorov于1933年提出的,它是针对单个总体的检验,检验样本是否来自指定分布f0.h0:f(x)=f0(x), h1:f(x)≠f0(x) (双侧检验)（1.63）
h0:f(x)≥f0(x), h1:f(x)h0:f(x)≤f0(x), h1:f(x)>f0(x) (单侧检验)（1.65）smirnov(斯米尔诺夫)拟合优度检验是smirnov于1939年提出的,它是针对两个总体的检验,检验两个总体的分布是否相同.h0:f(x)=g(x), h1:f(x)≠g(x) (双侧检验)（1.66）
h0:f(x)≥g(x), h1:f(x)h0:f(x)≤g(x), h1:f(x)>g(x) (单侧检验)（1.68）可以将smirnov检验看成kolmogorov检验的两样本情形,因此,将两个检验统称为kolmogorov-smirnov检验,简称k-s检验.
在r中,用ks.test()函数作k-s检验,其使用格式为ks.test(x, y, ...,
alternative = c("two.sided", "less", "greater"),
exact = null)参数x为数值向量.y或者是数值向量(smirnov检验),或者是描述分布函数的字符串(kolmogorov检验)....描述特定分布参数的字符串.
alternative为备择假设选项,取"two.sided"(默认值)表示双侧检验;取"less"表示备择假设为“<”的单侧检验;取"greater"表示备择假设为“>”的单侧检验.
exact为逻辑变量(默认值为null),表示是否精确计算p值.
例1.24 对一台设备进行寿命检验,记录10次无故障工作时间,并按从小到大的次序排列如下(单位：小时)：
420 500 920 1380 1510 1650 1760 2100 2300 2350
试用kolmogorov-smirnov检验方法检验此设备的无故障工作时间是否服从指数分布,且平均无故障工作时间为1500h.
解 指数分布的参数为λ,均值为1/λ,因此,平均时间为1500h,则参数为λ=1/1500.输入数据,调用ks.test()函数计算(程序名:exam0124.r)，p值(=0.2654)>0.05,无法拒绝原假设,因此,认为此设备的无故障工作时间服从指数分布,且平均无故障工作时间为1500h.> x <- c( 420, 500, 920, 1380, 1510, 1650, 17<code>`</code>javascript
60, 2100, 2300, 2350)
ks.test(x, "pexp", 1/1500)
one-sample kolmogorov-smirnov test
d = 0.3015, p-value = 0.2654
can("exam0125.data", nlines = 3)
y <- scan("exam0125.data", skip = 3)
ks.test(x, y)
two-sample kolmogorov-smirnov test
data: x and y
d = 0.23, p-value = 0.5286
alternative hypothesis: two-sidedp值(=0.5286)>0.05,接受原假设,即认为两个总体的分布函数是相同的.
kolmogorov-smirnov检验与pearson χ2检验相比,不需将样本分组,少了一个任意性,这是其优点.其缺点是只能用在理论分布为一维连续分布且分布完全已知的情形,适用面比pearson χ2检验小.研究也显示:在kolmogorov-smirnov检验可用的场合下,其功效一般来说略优于pearson χ2检验.
1.4.3 正态性检验
正态性检验是非参数检验中最常用的检验,除了用前面介绍的两种方法作检验外,还用许多检验方法,这里介绍shapiro-wilk(夏皮罗威尔克)正态性检验,它是利用来自总体x的样本x1,x2,…,xn,检验h0:总体x具有正态分布由于在检验中用w统计量作正态性检验,因此这种检验方法也称为正态w检验方法.
在r中,shapiro.test()函数完成shapiro-wilk正态性检验,其使用格式为shapiro.test(x)参数x是由样本构成的向量,并且向量的长度为3～5000.
例1.26 用shapiro-wilk正态性检验方法检验例1.1中考试成绩是否服从正态分布。
解 读取数据,调用shapiro.test()函数计算(程序名:exam0126.r).x <- scan("exam0101.data"); shapiro.test(x)
shapiro-wilk normality test
w = 0.8633, p-value = 0.0009852p值(=0.0009852)0.05,拒绝原假设,即考试成绩不服从正态分布,这与前面的分析(见例1.3和例1.4)是相符的.