非参数统计(R软件)参考答案讲述讲解

A.3, A.10, A.12

A.3 上机实践：将MASS数据包用命令library(MASS)加载到R中，调用自带“老忠实”喷泉数据集geyer，它有两个变量：等待时间waiting和喷涌时间duration，其中… (1) 将等待时间70min以下的数据挑选出来;

(2) 将等待时间70min以下,且等待时间不等于57min的数据挑选出来; (3) 将等待时间70min以下喷泉的喷涌时间挑选出来; (4) 将喷涌时间大于70min喷泉的等待时间挑选出来。 解:读取数据的R命令：

library(MASS);#加载MASS包 data(geyser);#加载数据集geyser

attach(geyser);#将数据集geyser的变量置为内存变量 (1) 依题意编定R程序如下：

sub1geyser=geyser[which(waiting<70),1];

#提取满足条件（waiting<70）的数据,which()，读取下标 sub1geyser[1:5];#显示子数据集sub1geyser的前5行 [1] 57 60 56 50 54

(2) 依题意编定R程序如下：

Sub2geyser=geyser[which((waiting<70)&(waiting!=57)),1]; #提取满足条件（waiting<70& (waiting!=57)的数据. Sub2geyser[1:5];#显示子数据集sub1geyser的前5行 [1] 60 56 50 54 60 ……

原数据集的第1列为waiting喷涌时间，所以用[which(waiting<70),2] (3)

Sub3geyser=geyser[which(waiting<70),2];

#提取满足条件（waiting<70）的数据,which()，读取下标 Sub3geyser[1:5];#显示子数据集sub1geyser的前5行 [1] 4.000000 4.383333 4.833333 5.450000 4.866667……

原数据集的第2列为喷涌时间，所以用[which(waiting<70),2] (4)

Sub4geyser=geyser[which(waiting>70),1];

#提取满足条件（waiting<70）的数据,which()，读取下标 Sub4geyser[1:5];#显示子数据集sub1geyser的前5行 [1] 80 71 80 75 77……. A.10

如光盘文件student.txt中的数据，一个班有30名学生，每名学生有5门课程的成绩，编写函数实现下述要求：

(1) 以data.frame的格式保存上述数据；

(2) 计算每个学生各科平均分，并将该数据加入(1)数据集的最后一列； (3) 找出各科平均分的最高分所对应的学生和他所修课程的成绩；

(4) 找出至少两门课程不及格的学生，输出他们的全部成绩和平均成绩；

(5) 比较具有(4)特点学生的各科平均分与其余学生平均分之间是否存在差异。 先将数据集读入R系统

student=read.table(\"…\

class(student):#显示数据集student的类型，

[1] \"data.frame\"#student是数据框

names(student);#显示数据框student的变量

[1] \"name\" \"math\" \"physics\" \"chem\" \"literat\" \"english\" \"mean\" #输出显示，数据框student有7个变量,第7个变量是平均值mean。 (1)

write.table(student,\"F:\\\\gzmu非参数统计\\\\data2014\\\\各章数据\\\\附录A\\\\x.txt\打开x.txt

\"name\" \"math\" \"physics\" \"chem\" \"literat\" \"english\" \"1\" \"Katty\" 65 61 72 84 79 \"2\" \"Leo\" 77 77 76 64 55 ……

(2) 依题意，要为原始数据集添加一个变量，即添加一列在最后。?[,6]=? me=rep(0,30);

for(i in 1:30){x=as.numeric(student[i,2:6]); me[i]=mean(x);} student$mean=me;

#上面程序的最后一行也可以如此:student[,7]=me names(student);

[1] \"name\" \"math\" \"physics\" \"chem\" \"literat\" \"english\" \"mean\" #如上显示，程序运行后数据框student添加了第7列mean.

(3) 依题意，在(2)的程序运行后做，要用到which(mean==max(mean))，如同A.3。 attach(student);

maxme=student[which(mean==max(mean)),];#找出最高平均分的记录，并赋予maxme; maxme;

name math physics chem literat english mean

15 Liggle 78 96 81 80 76 82.2

(4) 依题意，要用到二重的for和if. 由原数据框geyser给data1赋值时要用到数据转换： #x=as.numeric(student[i,2:6]);#读取student第i行2:6列的数据，#data1[k,]=x;#将x赋给data4 #的第k行。sum(x<60)是不及格门数。 Data1=student[1,];#赋初值 k=0;

for(i in 1:30){x=as.numeric(student[i,2:6]);

if (sum(x<60)>1){k=k+1;data1[k,]=student[i,];}} data1

name math physics chem literat english mean

1 Ricky 67 63 49 65 57 60.2 7 Simon 66 71 67 52 57 62.6 9 Jed 83 100 79 41 50 70.6 10 Jack 86 94 97 51 55 76.6 12 Jetty 67 84 53 58 56 63.6 13 Corner 81 62 69 56 52 64.0 14 Osten 71 64 94 52 52 66.6 25 Amon 74 79 95 59 59 73.2

(5) 依题意，要创造两个子集data4和data2, 用两样本的比较方法比较他们的平均成绩是否有显著差异。类似创造data1的方法，创造data2。并设x=data1$mean,y=data2$mean,比较二样本x,y是否有显著差异，由于还没有学非参数检验，试用t检验检验之(R的t检验函数为t.test(x,y)，原假设H0是两样本的均值相等，备择假设H1是两样本不等)。如果P值p-value<0.05,则拒绝原假设。 data2=student[1,];k=0;

for(i in 1:30){x=as.numeric(student[i,2:6]);

if (sum(x<60)<2){k=k+1;data2[k,]=student[i,];} };

下面做t检验

x=data1$mean;y=data2$mean; t.test(x,y)

Welch Two Sample t-test data: x and y

t = -3.0236, df = 9.309, p-value = 0.01386

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -11.493236 -1.684037

sample estimates: mean of x mean of y ： 67.17500 73.76364

结论：p-value = 0.01386<0.05,拒绝原假设，即认为两样本的平均成绩有显著差异。

A.12 在一张图上，用取值（-10，10）之间间隔均等的1000个点，采用不同的线型一颜色给制sin(),cos(),sin()+cos()的函数图形，图形要求有主标题和副标题，标示出从坐标 x=seq(-10,10,length=50);#构造向量x, x[1:5];#显示x的前5个数据

[1] -10.00000 -9.97998 -9.95996 -9.93994 -9.91992 sin=sin(x);#计算sin函数值 cos=cos(x);

sc=sin(x)+cos(x);

plot(sin~x,xlab=\"x\lines(cos~x,type=\"b\点线图 lines(sc~x,type=\"oitle(\"三角函数图\");

所得图形如下图，sin为黑色，cos为红色，sin+cos为绿色：

内容：

1.1; 1.2; 1.11;(附加题：1.4; 1.7; 1.8 有能力的可做附加题)

1.1 某批发市场从厂家购置一批灯泡，根据合同的规定，灯泡的使用的寿命平均不低于1000h。已知灯泡的使用寿命服从正态分布，标准差是20h，从总体中随机抽取了100只灯泡，得知样本均值为996h，问题是：批发商是否应该购买该批灯泡？

(1) 零假设和备择假设应该如何设置？给出你的理由。

(2)在零假设1000之下，给出检验的过程并做出决策，如果不能拒绝零假设，可能是哪里出了问题。 解：(1) 根据题意，问题的假设为H0:1000H1:1000

理由：1000是批发商的意愿，违背这个意愿，也就是拒绝原假设H0，他就购这批灯泡了。不能轻易否定的事情应置于被保护地位H0。这个问题的检验统计量为

ZX1000,z=(996-1000)/2=-2

20/nP值pvalue=pnorm(z,0,1)= 0.02275013, 在alpha=0.05时拒绝原假设，根据合同，不购这批灯泡。 (2) 假设检验问题：H0:0H1:0。

这样的假设是有问题的。假设检验是一种这样哲学：不轻易否定旧过程，置旧过程为H0于被保护的位置，而以小概率否定之。而一但被拒绝，以小概率事件原理，拒绝域不是小概率。反证H0不真。所谓“天欲报之，必先厚之”也，以显我为人之厚道，虽如此也不能保护H0,怪不得我也。面此假设违返旧过程，这样的假设毫无意义。

如果按照这个检验问题，检验的P值是pvalue=1- pnorm(z,0,1)= 0.9772499, 没有充分的理由拒绝原假设，结论也是不购进这批灯泡。但是犯批II类错误的概率是多少，鬼才知道呢。

1.2 考虑下面检验问题(不用计算已给的数据).

(1)如果X服从N(0, 1)分布，假设检验问题H0:0H1:1000。可以知道0.05的似然比

检验，如果X>1.645, 则将会拒绝H0： ，而且按照Neyman-Pearson引理，该检验是最优的。

现在，如果我们观察到X=2.1，该水平0.05的最优检验告诉我们拒绝=0的零假设，接受=1000的备择假设，你觉得有问题吗？问题在哪里？如何解决？ 答：有问题。假设检验在原假设条件成立下，得到拒绝域受

。而

只是其中的一种情况，故不能接受

，意思是拒绝。

，接

改进方法：可直接提出假设“均值为1000”进行检验。即检验

(2) 有两组学生的成绩，第一组为11名，成绩为x:100,99,99,100,100,100,100,99, 100, 99, 99; 第二组为2名，成绩为y: 50, 0. 我们对这两组数据作同样水平= 0.05的t检验（假设总体的均值为），H0:100H1:100。

对第二组数据的检验结果为：df=10, t= -2.8868,mean(x)= 99.54545, 单边检验(<100, less)的P值为p-value = 0.008099。所以拒绝原假设，认为<100。

对第二组数据检验的结果为：df=1, t值为-3，单边(<100, less)的P值为p-value = 0.1024，不拒绝原假设=100。但是mean(y)=25.

解：两个结论都不是合理的，t检验是针对正态数据做的，第一组数据事实上是两点分布，x的取值域为{99，100}，所以t检验的基本假设不满足，所以第一个检验是不合理的；第二组数据的t检验也是不合理的，样本量太少，不具有代表性。

(3)写出上面所用的t检验统计量，及p值的定义，解释水平=0.05的意义（注意，这里是一般情况，不要联系(2)中的具体数据例子），如果没有给定水平，如何用p值来做出结论？ 解：设样本X1,X2,...,Xn iid N(,2), 对于三种假设（双边假设，两个单边假设）都用同一个t统计量tX，p值p_value=P|T|t(双边检验，alternative=”two.side”)，S/np_value=PTt(右边检验, alternative=”greater”)，p_value=PTt(左边检验alternative=”less”)，其中T~t(n1)。p_value小于检验水平时拒绝原假设，接受H1 。则有 I. 双边假设检验H0:0H1:0，拒绝原假设H0tt/2 p_value=P|T|t< II. 右尾假设检验H0:0H1:0,拒绝原假设H0tt  p_value=PTt< III. 左尾假设检验H0:0H1:0,拒绝原假设H0tt  p_value=PTt< (4) 写出和t检验有关的关于均值的100(1-)%置信区间（不要联系(2)中的数据，说明你所有的符号的意义(如果有的话)）

解：t检验是在正态样本条件下做。确实，双边假设的t检验与置信区间一一对应。其双边假设检验式，有P|T|t≤P|T|t1

SSPXt/2Xt/21

nn其中随机变量T服从t(n-1)分布。S是正态样本的样本方差。

(5) 如果X1,X2,...X,n服从正态分布N(,2)，其中未知，写出有关的关于均值的100(1-)%的置信区间。一般来说，如果知道X1,X2,...,Xn有未知均值和已知方差2，但分布不知道，我们不能用上面写的置信区间？如果能，需要什么条件？根据是什么？用公式说明。 解：①如果X1,X2,...,Xn服从正态分布N(,2)，其中未知，写出有关的关于均值的100(1-)%的置信区间。用到下面两个统计量：

ZXX~N(0,1)，t~t(n1) /nS/n如果方差2已知，则用正态置信区间，用Z构造置信区间。

PXz/2Xz/21

nn②如果方差2未知，则用t构造置信区间：

SSPXt/2Xt/21

nn③如果知道X1,X2,...,Xn有未知均值和已知方差2，但分布不知道，我们不能用上面写的置信区间，用切比雪夫不等式构造置信区间：

22P|X|12，令=2

nn2PXX12

n2(6)在切比雪夫不等式中，令B=2，12，，所以对给定的检验水平，

nn1B2

n1

1.11 (数据光盘文件：beenswax.txt)为探测蜂蜡结构，生物学家做了很多实验，在每个蜂蜡里碳氢化合物(hydrocarbon)所占的比例对蜂蜡结构有特殊的意义，数据中给出了一些观测。 (1)画出beenswax数据的经验累积分布、直方图和Q-Q图。 (2)找出0.9,0.75,0.50,0.25,0.10的分位数。 (3)这个分布是高斯分布吗？ 解：

beenswax=read.table(\"F:\\\\gzmu非参数统计\\\\data2014\\\\各章数据\\\\第1章\\\\beenswax.txt\attach(beenswax); names(beenswax)

[1] \"MeltingPoint\" \"Hydrocarbon\"

说明beenswax有两个变量：\"MeltingPoint\" \"Hydrocarbon\"，分别表示，熔点和碳氢化合物所占比例。

(1) 依题意，对Hydrocarbon的作图程序如下得图1.11-1 cdf=ecdf(Hydrocarbon);#计算经验分布函数 par(mfrow=c(2,2));#定义图矩阵为2行2列 plot(cdf);hist(Hydrocarbon);

qqnorm(Hydrocarbon);qqline(Hydrocarbon)

图1.11-1 图1.11-2

将上述程序中的Hydrocarbon替换成MeltingPoint,对MeltingPoint的作图程序如下得图1.11-2 cdf=ecdf(MeltingPoint);#计算经验分布函数 par(mfrow=c(2,2));#定义图矩阵为2行2列 plot(cdf);hist(MeltingPoint);

qqnorm(MeltingPoint);qqline(MeltingPoint)

(3) 从直方图看，两者基本成对称，钟形，从两者的正态Q-Q图，也知道，两者的散点基本在两条直线的附近。所以两近似正态分布（高斯分布）。

对Hydrocarbon和MeltingPoint做ks.test,P值分别为：0.9766和0.7774, 两个检验都没有拒绝原假设（数据呈正态分布）。程序如下：

ks.test(Hydrocarbon,pnorm,mean(Hydrocarbon),sd(Hydrocarbon)); ks.test(MeltingPoint,pnorm,mean(MeltingPoint),sd(MeltingPoint));

内容：

2.1, 2.2, 2.4, 2.1 2, 2.14

2.1 超市经理想了解每位顾客在该超市购买的商品平均件数是否为10件，随机观察12位顾客，得到如下数据： 顾客 1 件数 22 2 9 3 4 4 5 5 1 6 16 7 15 8 26 9 47 10 8 11 31 12 7 (1)采用符号检验进行决策。

(2)采用Wilcoxon符号秩检验进行决策，比较它和符号检验的结果。 （如果分布对称，则Wilcoxont符号秩检验较优，P值小者较优） 解：(1) 采用符号检验进行决策： 根据题意，检验的假设为双边假设

H0:median10H1:median10

x=c(22,9,4,5,1,16,15,26,47,8,31,7);

sg=sum(x>10);sl=sum(x<10);n1=sg+sl;k=min(sg,sl); binom.test(k,n1,0.5); 结果输出：

Exact binomial test data: k and n1

图2.1.1 数据分布直方图 number of successes = 6, number of trials = 12, p-value = 1

alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.2109446 0.7890554 sample estimates:probability of success 0.5 p-value = 1，不拒绝原假设H0

(2) Wilcoxon符号秩检验，假设如果(1)：

Wilcoxon signed rank test with continuity correction data: x - 10

V = 53, p-value = 0.2892

alternative hypothesis: true location is not equal to 0 p-value = 0.2892, 没有充分理由拒绝原假设。

注：虽然两个检验的结论相同，但我们认为(1)可靠。因为数据的分布不是对称，而后者是基于对称分布的。而本题的数据分布直方图如下，显然是不对称的，所针对本题数据，wilcox.test不可靠。

2.2考查某疾病的患者共计350名，男性150人，女性200人，问该疾病得病的男女性别比是否为1:1，即其男女比例是否各为1/2？

提示：用中心极限定理，正态近似检验，即Demoive-Laplace中心极限定理:p=0.5,n=350, X~b(350,0.5),E(X)=175, Var(X)=npq=n/4=350/4。标准化X近似于标准正态。 解：根据题意，设男性患者的比例为p,则检验的假设为H0:p0.5H1:p0.5

设男性患者数为X，则X~b(350,0.5),E(X)=175, Var(X)=npq=n/4=350/4。标准化X近似于标准正态。

ZXnp150175-2.672612 ~N0,1，zn/4npqp-value=2*min(pnorm(z,0,1),1- pnorm(z,0,1))= 0.007526315, 拒绝原假设p=0.5，认为患者中男性比率不是0.5,

男女比例不是1:1.

注：究其实，男性患者的比率显著地<0.5.

2.4 下表中的数据是两个篮球联赛中三分球的进球次数，该数据的目的是考察两个联赛三分球得分次数是否存在显著性差异。 (1)符号检验。

(2)配对Wilcoxon符号秩检验。

(3)在这些数据中哪个检验更好？为什么?（P值小者好） 三分球进球次数 队伍序号 1 2 3 4 5 6 7 8 9 10

联赛1 91 46 108 99 110 105 191 57 34 81 联赛2 81 51 63 51 46 45 66 64 90 28 解：设联赛1和联赛2的三分球得分次数分别为X和Y，题意只问“X和Y”是否存在显著差异，所以检验的假设为H0:meXmeYH1:meXmeY

设Z=X-Y，问题转化为H0:meZ0H1:meZ0

(1) 检验的R程序为：

x=c(91,46,108,99,110,105,191,57,34,81);y=c(81,51,63,51,46,45,66,64,90,28); z=x-y;

sg=sum(z>0);sl=sum(z<0);n1=sg+sl;k=min(sg,sl) binom.test(k,n1,0.5)

Exact binomial test data: k and n1

number of successes = 3, number of trials = 10, p-value = 0.3438 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.06673951 0.65245285 sample estimates:probability of success 0.3

P值p-value = 0.3438,不拒绝原假设，认为两个联赛的三分球得分次数没有显著差异。 (2) 作z的直方图如图2.4.1，图形显示z的分布不存在显著不对称的迹象，可以做wilcox.test wilcox.test(z)

Wilcoxon signed rank test data: z

V = 45, p-value = 0.08398

alternative hypothesis: true location is not equal to 0

检验的P值p-value = 0.08398,在alpha=0.05下，不拒绝原假设。与符号检验的结

图2.4.1 z的直方图 论相同，但P值小了很多。

(3) 在如上的检验中，由于数据的分布不存在显著不对称的迹象，wilcox.test是可靠的，因而wilcox.test理好。事实wilcox.test的P值小了很多，更能区分差异。在检验可靠的情形下，P值越小越好。

2.12 在白令海所捕捉的12岁的某种鱼的长度(单位：cm)样本为 长度/cm 64 65 66 67 68 69 70 71 72 73 74 75 77 78 79 数目 1 2 1 1 4 3 4 5 3 3 0 1 6 1 1 您能否同意所声称的12岁的这种鱼的长度的中位数总是在69~72cm之间？ 解：这是求置信区间的问题，设=0.05.

x=c(64,65,65,66,67,68,68,68,68,69,69,69,70,70,70,70,71,71,71,71,71,72,72,72,73,73,73,75,77,77,77,77,77,77,78,83);

数据探索：正态Q-Q图和密度函数图如下

两者显示数据x近似于对称分布，ks正态性检验的P值为0.58，也没有拒绝正态性假设，因此可以认为数据分布不拒绝对称性假设。因此可以做Walsh中位数置信区间，基于Bootstrap方差估计的中位数正态置信区间、枢

轴量置信区间、分位数置区间，下面求walsh置信区间。 (1) walsh中位数置信区间 walsh=NULL;n=length(x);

for(i in 1:n){for (j in i:n){w=(x[i]+x[j])/2; walsh=c(walsh,w);}}

list(med=median(walsh), nwalsh=length(walsh));

# median(walsh)=71, length(walsh)=666 #编程求walsh中位数的(1-)*100%=95%的置信区间 walsh.conf=function(x,alpha){ walsh=NULL;n=length(x);

for(i in 1:n){for (j in i:n){w=(x[i]+x[j])/2; walsh=c(walsh,w);}}

nw=length(walsh); #walsh的长度

walsh.sort=sort(walsh);#搜索walsh中位数的置信区间，对称地砍掉左尾和右尾 for(k in seq(1,(nw/2),1)){F=pbinom(nw-k,nw,0.5)- pbinom(k,nw,0.5); if (F<(1-alpha)){lk=k-1;break}}

lci=walsh.sort[lk];uci=walsh.sort[nw-lk+1];

list(lci=lci,uci=uci,lk=lk,uk=nw-lk)} #调用函数 walsh.conf(x,0.05) $lci= 71, $uci=71.5

结论：12岁的这种鱼的长度的中位数的95%的walsh置信区间是（71, 71.5）(cm). (2) 其它置信区间，基于Bootstrap方差的枢轴区间是最好的，它是（69，73），还是没有Walsh区间好，因为数据分布是对称的。

依walsh平均，可以说12岁的这种鱼的长度在69~72之间（置信水平≥95%）。

2.14 社会学家欲了解抑郁症的发病率是否在一年时间随季节的不而不同，他使用了来年一所大医院的病人数据，按一个4个季节,依次记录过去5年中第一次被确诊为患抑郁症的病人数，数据如下表(单位：人)

季节 春季 夏季 秋季 冬季 合计 人数 495 503 491 581 2070 请问：发病率是否与季节有关？

解：这是一个假设问题。也称为独立性检验问题。如果两者独立，即无关，则发病人数在4个季节是均匀（发病率为1/4），否则两者是相关的。Pearson检验过程如下：

H0；p1=p2=p3=p4=1/4；H1；p1，p2，p3，p4不全等； V=c(495,503,491,581);p=1/4;n=sum(V);df=4-1; chi2=sum((V-n*p)^2/(n*p))

pvalue=1-pchisq(chi2,df);pvalue;#请思考：为什么用右尾概率？ [1]0.01453647

结论：在=0.05时拒绝原假设，认为发病率与季节有关。具体地说，冬天的发病率高（p3= 0.2807）。当然，为了要得到科学的结论，应该要规范抽样，使得样本有代表性，毕竟一个医院的数据其代表性是值得商榷的。

内容

P106: 3.1; 3.4; 3.5.

3.1 在一项研究毒品对增强人体攻击性影响的实验中，组A使用安慰剂，组B使用毒品，试验后进行攻击性测试，测量得分显示在如下表中(得分越高表示攻击性越强) 组A 10,8,12,16,5,9,7,11,6 2组B 12,15,20,18,13,14,9,16 (1)给出这个实验的零假设. (2)画出表现这些数据的曲线图.

(3)分析这些数据用哪种检验方法最合适. (4)用您选择的检验对数据进行分析.

(5)是否有足够的证据拒绝零假设？如何解释数据？

解：(1)这个实验的目的是要检验毒品是否具有显著的攻击性。根据假设检验的原则，其零假设其位置参数（均值或中位数）是无显著差异，即检验假设为：.

H0:MAMBH1:MAMB

(2)

A=c(10,8,12,16,5,9,7,11,6); B=c(12,15,20,18,13,14,9,16);

min=min(c(A,B));max=max(c(A,B));

plot(A,type=\"b\lines(B,type=\"bitle(\"数据A,B折线图\"); 折线图如图3.1.1.

能更好地反映数据还有箱线图，程序如下，图如图3.1.2 group=factor(rep(c(\"A\plot(c(A,B)~group)

图3.1.1 数据A、B折线性图 图3.1.2 数据A、B箱线图

从图看，药品B的攻击性是乎强一些，有否显著地强，有待于检验。

(3)如果两样本都呈正态分布，可以进行二样本t检验，如果两样本分布相似，可进行Wilcoxon秩和检验。 二样本正态性检验的程序和结果如下： ks.test(A,pnorm,mean(A),sd(A))

One-sample Kolmogorov-Smirnov test data: A

D = 0.1047, p-value = 0.9997 alternative hypothesis: two-sided

因为检验的P值为0.9997,没有充分的理由拒绝A的正态性假设。 ks.test(B,pnorm,mean(B),sd(B))

One-sample Kolmogorov-Smirnov test data: B

D = 0.0991, p-value = 1

alternative hypothesis: two-sided.

因为检验的P值为1,没有充分的理由拒绝B的正态性假设。

所以可以进行t检验

t.test(A,B,alternative=\"less\ Welch Two Sample t-test data: A and B

t = -3.1763, df = 14.686, p-value = 0.0032

alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf -2.366982

sample estimates:mean of x mean of y 9.333333 14.625000 再做两样本分布相似检验

ks.test(A-median(A),B-median(B))

Two-sample Kolmogorov-Smirnov test data: A - median(A) and B - median(B) D = 0.1389, p-value = 0.9998 alternative hypothesis: two-sided

因为检验的P值为0.9998,没有充分理由拒绝两样本分布相似的假设，所以可做wilcox.test wilcox.test(A,B,alternative=\"less\")

Wilcoxon rank sum test with continuity correction data: A and B

W = 9.5, p-value = 0.006097

alternative hypothesis: true location shift is less than 0

因为t检验的P值为0.0032,而Wilcoxon秩和检验的P值为0.00609,在=0.01时，两者均有充分的理由拒绝零假设，认为毒品B具有显著的攻击性。

(4) (5)因为t检验的P值为0.0032,而Wilcoxon秩和检验的P值为0.00609,在=0.01时，两者均有充分的理由拒绝零假设，认为毒品B具有显著的攻击性。

3.4 两个不同学院教师一年的课时量分别为（单位：学时）

A学院：321,266,256,386,330,329,303,334,299,221,365,250,258,342,243,298,238,317 B学院：488,593,507,428,807,342,512,350,672,589,665,549,451,492,514,391,366,469

根据这两个样本，两个学院教师讲课的课时是否存在显著差异？估计这些差异。从两个学院教师讲课的课时来看，教师完成讲课任务的情况是否类似？给出检验和判断。 提示：先检验“教师完成讲课任务的情况是否类似”，再选择检验方法，推断是否存在显著差异。 解：

A=c(321,266,256,386,330,329,303,334,299,221,365,250,258,342,243,298,238,317); B=c(488,593,507,428,807,342,512,350,672,589,665,549,451,492,514,391,366,469); (1) 检验“教师完成讲课任务的情况是否类似” 方法ks.test检验：

H0:两样本分布类似H1:两样分布不类似。

ks.test(A-median(A),B-median(B)) 检验结果：

Two-sample Kolmogorov-Smirnov test data: A - median(A) and B - median(B) D = 0.2778, p-value = 0.5026 alternative hypothesis: two-sided

因为检验的P值为0.5026,不拒绝零假设H0，即不拒绝两样本分布类似的假设。

注：如果分别用正态性检验，则在不拒绝正态性假设的基础上，还要检验两样本方差齐性。思考一下为什么？

(2)在(1)的检验中，两样本分布相似，所以可以用Wilcoxon秩和检验检验两样本中位数是否有显著差异：

H0:MedxMedyH1:MedxMedy

wilcox.test(A,B) 检验结果：

Wilcoxon rank sum test with continuity correction data: A and B

W = 5.5, p-value = 7.977e-07

alternative hypothesis: true location shift is not equal to 0

因为检验的P值为7.977e-07<<0.01,所以拒绝零假设，两样本的中位数有显著差异。两学院教师的教学任务有显著差异。

(3)可以在(2)的基础上进一步检验，两样本A与B不但分布相似，而且相似于正态分布（两者均呈正态分布），所以可以用二样本t检验： t.test(A,B) 检验结果：

Welch Two Sample t-test data: A and B

t = -6.8841, df = 21.916, p-value = 6.637e-07

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -276.8201 -148.6243

sample estimates: mean of x mean of y 297.5556 510.2778

因为检验的P值为6.637e-07<<0.01,所以拒绝零假设，两样本的均值存在显著差异。即两学院教师的平均教学课时有显著差异。两学院教师平均教学课时分别为297.5556，510.2778，B学院教师的平均课时显著地高于A学院。

(4)两样本位置(均值、中位数)差的各种估计，置信区间讨论：…… (5)两样本密度估计，非参数密度估计：……

您有什么想法？将这些方法实施于理学院本科学生成绩分析，教师教学任务的统计分析？您愿意做这些平凡的实际工作？如果您展开充分的思考，提升到社会学乃至心理学，我看是可以做学位论文的。精彩的统计分析工作还可以在以后章节遇到。世界上怕就怕您高不成低不就啊！

3.5 对A和B两块土壤有机质含量抽检结果如下，试用Mood和Moses两种方法检验两组数据的方差是否存在显著差异。 A B 8.8,8.2 5.6,4.9 8.9,4.2 3.6,7.1 5.5,8.6 6.3,3.9 13.0,14.5 16.5,22.6, 20.7,19.6, 18.4,21.3, 24.2,19.6, 11.7, 18.9,14.6, 19.8,14.5 解：

A=c(8.8,8.2,5.6,4.9,8.9,4.2,3.6,7.1,5.5,8.6,6.3,3.9);

B=c(13.0,14.5,16.5,22.6,20.7,19.6,18.4,21.3, 24.2,19.6,11.7,18.9,14.6,19.8,14.5);

(1) Mood方差检验是数据中心化后，用混合样本的秩代替离差平方和公式中的原始数据。即设样本

X(X1,X2,...,Xm)iidF(x/1)，Y(Y1,Y2,...,Yn)iidF(x/2)，检验的假设为

22。 H0:1222H1:122再设X在混合样本c(X,Y)中的秩为R=(R1,R2,…,Rm)，当H0成立时，混合样本c(X,Y)= (X1,X2,...,Xm;

Y1,Y2,...,Yn) iid Fx/

ERi而秩统计量

mnimn1 2i1mnMRi(mn1)/2

2i1m应该不大(Ri在平均值(m+n+1)/2附近波动),而当X的方差大于Y的方差时，Ri会在远离平均值(m+n+1)/2的地方出现，因而当M超大时，拒绝零假设。

检验可以编程计算，也可以调用R的现成函数mood.test()。本题数据运行 mood.test(A,B)

Mood two-sample test of scale data: A and B

Z = 0.6342, p-value = 0.526 alternative hypothesis: two.sided

由于P值为0.526，没有充分理由拒绝原假设。

(2) Moses的方法是将两样本分组，用各组的离差平方和反映方差。分组要注意到每组中至少有3个样本。本题中样本容量分别为12，15，所以分别分为4组，5组。 SSA=NULL;

for (i in 1:4){group=A[((i-1)*3+1):(3*i)];SSA=c(SSA,2*var(group))} SSA=(5.786667 12.860000 6.140000 11.046667) SSB=NULL;

for (j in 1:5){group=B[((i-1)*3+1):(3*i)];SSB=c(SSB,2*var(group))} SSB=(38.24667 38.24667 38.24667 38.24667 38.24667); wilcox.test(SSA,SSB)

Wilcoxon rank sum test with continuity correction data: SSA and SSB W = 0, p-value = 0.0108

alternative hypothesis: true location shift is not equal to 0

结论：两组数据的方差有显著差异，由median(SSA)=8.593333,median(SSB)=38.24667,所B数据方差显著大于A组数据的方差。

Moses的缺点是，分组后样本量缩小了，很不好。用Bootstrap方法，直接比较Bootstrap样本的方差，思想方法简单：重抽样B次，各得方差的B个大样本，由大样本理论比较两样本方差。 (3) Bootstrap方法

x=c(8.8,8.2,5.6,4.9,8.9,4.2,3.6,7.1,5.5,8.6,6.3,3.9);

y=c(13.0,14.5,16.5,22.6,20.7,19.6,18.4,21.3, 24.2,19.6,11.7,18.9,14.6,19.8,14.5); VBx=NULL;VBy=NULL;

nx=length(x);ny=length(y);B=1000;

for (i in 1:B){xb=sample(x,nx,T);Vbx=var(xb); VBx=c(VBx,Vbx);

yb=sample(y,ny,T);Vby=var(yb);

VBy=c(VBy,Vby);};

MVx=mean(VBx);MVy=mean(VBy); Varxy=var(VBx)+var(VBy);

Z=(MVx-MVy)/sqrt(Varxy);#计算Z值，大样本就是要用Z值，要用中心极限定理。 p1=pnorm(Z);p2=1-pnorm(Z); pvalue=2*min(p1,p2);pvalue; [1] 0.01107871

此结果与（2）相同。

内容：P143 4.1；P144-4.4; P144-4.5;

4.1对A,B,C三个灯泡厂生产的灯泡进行寿命测试，每种品牌随机试验不等量灯泡，结果得到如下列寿命数据(单位：天)，试比较三品牌灯泡寿命是否相同。 A B C 83 64 67 62 70 85 81 80 78 88 89 79 90 95 解：（1）三个样本A、B、C均为独立随机样本，非区组试验数据，样本量不同，只能用Kruskal秩方差分析方法。检验的假设

H0:三样本的中位数相同 VS H1:三样本的中位数不全相同

因为样本少，免做样本数据分布相似检验，直接做kruskal.test,程序如下： A=c(83,64,67,62,70);B=c(85,81,80,78);C=c(88,89,79,90,95); n1=length(A); n2=length(B); n3=length(C); x=c(A,B,C);

group=factor(rep(1:3,c(n1,n2,n3))); kruskal.test(x~group) 结果：

Kruskal-Wallis rank sum test data: x by group

Kruskal-Wallis chi-squared = 7.8229, df = 2, p-value = 0.02001

即检验的P值为0.02001，拒绝原假设，即A,B,C三个灯泡厂生产的灯泡的寿命的中位数有显著差异。 （2）进一步分析差异出自何处，请看箱线盒须图： plot(x~group)

图4.1-1 三个厂的灯泡寿命

图4.1-1显示，至少处理C与处理A有显著差异，由于灯泡寿命是望大的，所以C厂生产的灯泡寿命最长,最优。

（3）两两比较的程序和结果如下

A=c(83,64,67,62,70);B=c(85,81,80,78);C=c(88,89,79,90,95); n1=length(A); n2=length(B); n3=length(C);k=3;

n=c(n1,n2,n3);alpha=0.05; alphas=alpha/(k*(k-1));Z=qnorm(alphas,0,1); N=sum(n);MST=N*(N+1)/12;

x=c(A,B,C);R=rank(x);Rbar=rep(0,k); group=factor(rep(1:3,c(n1,n2,n3)));

for(i in 1:k){Rbar[i]=median(R[group==i])} d=NULL;

for(i in 1:(k-1)){for (j in (i+1):k){SE=sqrt(MST*(1/n[i]+1/n[j])); d=c(d,abs(Rbar[i]-Rbar[j])/SE)}} nd=length(d);dsig=rep(0,nd)

for (i in 1:nd){if (d[i]<=Z){dsig[i]=1}}

dsig;#dsig=0,两者有显著差异。length(dsig)=k*(k-1)/2 dsig 1-2 1-3 2-3 [1] 0 0 0

说明多重比较中，两两均有显著差异。

4.4 下表是美国三大汽车公司(A,B,C三种处理)的五种不同的车型某年产品的油耗，试分析不同公司的油耗

是否存在差异。 A B C 1 20.3 25.6 24.0 2 21.2 24.7 23.1 3 18.2 19.3 20.6 4 18.6 19.3 19.8 5 18.5 20.7 21.4 解：（1）事实上，这张表的实验数据是双因素(公司，车型)试验数据表，A、B、C的样本数据独立但不同分布，所以要检验不同公司的车的油耗，即检验A、B、C的差异，要剔除区组之影响，不能用kruskal.test,

而只能用Friedman.test,检验程序如下：

A=c(20.3,21.2,18.2,18.6,18.5); B=c(25.6,24.7,19.3,19.3,20.7); C=c(24.0,23.1,20.6,19.8,21.4);

n1=length(A); n2=length(B); n3=length(C); x=c(A,B,C);

M=matrix(x,3,5,byrow=T);

friedman.test(t(M)) ;#这里要小心啊，该检验以区组为列 检验结果：

Friedman rank sum test data: t(M)

Friedman chi-squared = 7.6, df = 2, p-value = 0.02237

即检验的P值为0.02237，拒绝原假设，即A,B,C三个公司汽车的油耗有显著差异。 （2）进一步分析差异出自何处，请看箱线盒须图： plot(x~group)

图4.4-1 三个公司汽车的油耗

图4.4-1显示，至少处理C与处理A有显著差异，由于汽车油耗是望小的，所以公司A汽车油耗是最少的，是最优的。如果从只从油耗考虑，买汽车应该买A公司的汽车。 （3）两两比较的程序与结果如下 A=c(20.3,21.2,18.2,18.6,18.5); B=c(25.6,24.7,19.3,19.3,20.7); C=c(24.0,23.1,20.6,19.8,21.4);

n1=length(A); n2=length(B); n3=length(C);k=3;

n=c(n1,n2,n3);alpha=0.05; alphas=alpha/(k*(k-1));Z=qnorm(alphas,0,1); N=sum(n);MST=N*(N+1)/12;

x=c(A,B,C);R=rank(x);Rbar=rep(0,k); group=factor(rep(1:3,c(n1,n2,n3)));

for(i in 1:k){Rbar[i]=median(R[group==i])} d=NULL;

for(i in 1:(k-1)){for (j in (i+1):k){SE=sqrt(MST*(1/n[i]+1/n[j])); d=c(d,abs(Rbar[i]-Rbar[j])/SE)}} nd=length(d);dsig=rep(0,nd)

for (i in 1:nd){if (d[i]<=Z){dsig[i]=1}}

dsig;#dsig=0,两者有显著差异。length(dsig)=k*(k-1)/2 #dsig 1-2 1-3 2-3 [1] 0 0 0

说明多重比较中，两两均有显著差异。

4.5 在一项健康试验中，有三种生活方式，它们的减肥效果如下表 生活方式 一个朋后 减少的质量 (单位：500g) ni 1 3.7 3.7 3.0 3.9 2.7 5 2 7.3 5.2 5.3 5.7 6.5 5 3 9.0 4.9 7.1 8.7 4 人们想知道的是从这些数据能否得出它们的减肥效果(位置参数)是一样的，如果效果不等，试根据上面这些数据选择方法检验哪一各效果最好，哪一种最差。 解：（1）根据试验设计和数据，这不是区组试验，但也不是简单随机样本(因为存在人体这个“混杂”因素)，只能勉强用Kruskal秩方差分析。检验程序和结果如下：

x1=c(3.7,3.7,3.0,3.9,2.7);x2=c(7.3,5.2,5.3,5.7,6.5);x3=c(9.0,4.9,7.1,8.7); n1=length(x1); n2=length(x2); n3=length(x3); x=c(x1,x2,x3);

group=factor(rep(1:3,c(n1,n2,n3))); kruskal.test(x~group) 结果输出：

Kruskal-Wallis rank sum test data: x by group

Kruskal-Wallis chi-squared = 9.4322, df = 2, p-value = 0.00895

即检验的P值为0.00895，拒绝原假设，即三种生活方式的减肥效果有显著差异。 （2）进一步分析差异出自何处，请看箱线盒须图：

plot(x~group)

图4.5-1 三种生活方式的减肥效果

图4.5-1显示，至少第1种生活方式与第3种生活方式的中位数有显著差异，由于减肥效果是望大的，所以第3种减肥效果是最优的。

（3）两两比较的程序与结果如下

A=c(3.7,3.7,3.0,3.9,2.7);B=c(7.3,5.2,5.3,5.7,6.5);C=c(9.0,4.9,7.1,8.7); n1=length(A); n2=length(B); n3=length(C);k=3;

n=c(n1,n2,n3);alpha=0.05; alphas=alpha/(k*(k-1));Z=qnorm(alphas,0,1); N=sum(n);MST=N*(N+1)/12;

x=c(A,B,C);R=rank(x);Rbar=rep(0,k); group=factor(rep(1:3,c(n1,n2,n3)));

for(i in 1:k){Rbar[i]=median(R[group==i])} d=NULL;

for(i in 1:(k-1)){for (j in (i+1):k){SE=sqrt(MST*(1/n[i]+1/n[j])); d=c(d,abs(Rbar[i]-Rbar[j])/SE)}} nd=length(d);dsig=rep(0,nd)

for (i in 1:nd){if (d[i]<=Z){dsig[i]=1}}

dsig;#dsig=0,两者有显著差异。length(dsig)=k*(k-1)/2 #dsig 1-2 1-3 2-3 [1] 0 0 0

说明多重比较中，两两均有显著差异。