非参数统计(第二版)习题R程序.doc

P37.例2.1
build.price<-c(36,32,31,25,28,36,40,32,41,26,35,35,32,87,33,35);build.price
hist(build.price,freq=FALSE)#直方图
lines(density(build.price),col="red")#连线
#方法一：m<-mean(build.price);m#均值
D<-var(build.price)#方差
SD<-sd(build.price)#标准差S
t=(m-37)/(SD/sqrt(length(build.price)));t#t统计量计算检验统计量
t=
[1] -0.1412332
#方法二：t.test(build.price-37)#课本第38页
例2.2
binom.test(sum(build.price<37),length(build.price),0.5)#课本40页
例2.3
P<-2*(1-pnorm(1.96,0,1));P
[1] 0.04999579
P1<-2*(1-pnorm(0.7906,0,1));P1
[1] 0.4291774
> 例2.4
> p<-2*(pnorm(-1.96,0,1));p
[1] 0.04999579
>
> p1<-2*(pnorm(-0.9487,0,1));p1
[1] 0.3427732

例2.5（P45）
scores<-c(95,89,68,90,88,60,81,67,60,60,60,63,60,92,
60,88,88,87,60,73,60,97,91,60,83,87,81,90);length(scores)#输入向量求长度
ss<-c(scores-80);ss
t<-0
t1<-0
for(i in 1:length(ss)){
if (ss[i]<0) t<-t+1#求小于80的个数
else t1<-t1+1求大于80的个数
}
t;t1
> t;t1
[1] 13
[1] 15
binom.test(sum(scores<80),length(scores),0.75)
p-value = 0.001436<0.01
Cox-Staut趋势存在性检验P47
例2.6
year<-1971:2002;year
length(year)
rain<-c(206,223,235,264,229,217,188,204,182,230,223,
227,242,238,207,208,216,233,233,274,234,227,221,214,
226,228,235,237,243,240,231,210)
length(rain)
1)该地区前10年降雨量是否变化？
t1=0
for (i in 1:5){
if (rain[i]x(k1),16,0.5)#= 0.0002593994
y1<-(1+16)/(2^16);y1#=0.0002593994
plot(year,rain)
abline(v=(1971+2002)/2,col=2)
lines(year,rain)
anova(lm(rain~(year)))
随机游程检验（P50）
例2.8
client<-c("F","M","M","M","M","M","F","M",
"M","F","M","M","M","M","F","M","F",
"M","M","M","F","F","F","M","M","M");client
n<-length(client);n
n1<-sum(client=="M");n1
n0<-n-n1;n0
t1<-0
for (i in 1:(length(client)-1)){
if (client[i]==client[i+1]) t1<-t1
else t1<-t1+1
}
R<-t1+1;R#=12
#find rejection region（不写）
rl<-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0));rl
ru<-2*n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0));ru#=15.33476（课本为ru=17）
例2.9
shuju39<-data.frame(read.table
("SHUJU39.txt",header=TRUE));shuju39
attach(shuju39)
sum.a=0
sum.b=0
sum.c=0
for (i in 1:length(id)){
if (pinzhong[i]=="A") sum.a<-sum.a+chanliang[i]
else if (pinzhong[i]=="B") sum.b<-sum.b+chanliang[i]
else fuhao<-sum.c<-sum.c+chanliang[i]
}
sum.a;sum.b;sum.c
ma<-sum.a/4
mb<-sum.b/4
mc<-sum.c/4
ma;mb;mc
fuhao<-rep("a",12);fuhao
for (i in 1:length(id)){
if (pinzhong[i]=="A" & ((chanliang[i]-ma)>0)) fuhao[i]<-"+"
else if (pinzhong[i]=="B" & ((chanliang[i]-mb)>0)) fuhao[i]<-"+"
else if (pinzhong[i]=="C" & ((chanliang[i]-mc)>0)) fuhao[i]<-"+"
else fuhao[i]<-"-"
}
fuhao

#利用上题编程解决检验的随机性
n<-length(fuhao);n
n1<-sum(fuhao=="+");n1
n0<-n-n1;n0
t1<-0
for (i in 1:(length(fuhao)-1)){
if (fuhao[i]==fuhao[i+1]) t1<-t1
else t1<-t1+1
}
R<-t1+1;R
#find rejection region
rl<-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0));rl
ru<-2n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0));ru
例2.10（P52）library(quadprog)# 不存在叫‘quadprog’这个名字的程辑包
library(zoo)# 不存在叫‘zoo’这个名字的程辑包
library(tseries)# 不存在叫‘tseries’这个名字的程辑包
run1=factor(c(1,1,1,0,rep(1,7),0,1,1,0,0,rep(1,6),0,rep(1,4),
0,rep(1,5),rep(0,4),rep(1,13)))；run1
y=factor(run1)
runs.test(y)# 错误 没有"runs.test"这个函数

Wilcoxon符号秩检验
W+在零假设下的精确分布
#下面的函数dwilxonfun用来计算W+分布密度函数，
即P(W+=x)的一个参考程序！
dwilxonfun=function(N){
a=c(1,1) #when n=1 frequency of W+=1 or o
n=1
pp=NULL #distribute of all size from 2 to N
aa=NULL #frequency of all size from 2 to N
for (i in 2:N){
t=c(rep(0i),a)
a=c(a,rep(0,i))+t
p=a/(2^i) #density of Wilcox distribut when size=N
}
p
}
N=19 #sample size of expected distribution of W+
y<-dwilxonfun(N);y
#计算P(W+=x)中的x取值的R参考程序！！
dwilxonfun=function(N){
a=c(1,1) #when n=1 frequency of W+=1 or o
n=1
pp=NULL #distribute of all size from 2 to N
aa=NULL #frequency of all size from 2 to N
for (i in 2:N){
t=c(rep(0,i),a)
a=c(a,rep(0,i))+t
p=a/(2^i) #density of Wilcox distribut when size=N
}
a
}
N=19 #sample size of expected distribution of W+
y<-dwilxonfun(N);length(y)-1
hist(y,freq=FALSE)
lines(density(y),col="red")
例2.12（P59）
ceo<-c(310,350,370,377,389,400,415,425,440,295,
325,296,250,340,298,365,375,360,385);length(ceo)
#方法一
wilcox.test(ceo-320)
#方法二
ceo.num<-sum(ceo>320);ceo.num
n=length(ceo)
binom.test(ceo.num,n,0.5)

例2.13(P61)
a<-c(62,70,74,75,77,80,83,85,88)
walsh<-NULL
for (i in 1:(length(a)-1)){
for (j in (i+1):length(a)){
walsh<-c(walsh,(a[i]+a[j])/2)
}
}
walsh=c(walsh,a)
NW=length(walsh);NW
median(walsh)
2.5单组数据的位置参数置信区间估计(P61)
例2.14‘
stu<-c(82,53,70,73,103,71,69,
80,54,38,87,91,62,75,65,77);stu
alpha=0.05
rstu<-sort(stu);rstu
conff<-NULL;conff
n=length(stu);n
for(i in 1:(n-1)){
for (j in (i+1):n){
conf=pbinom(j,n,0.5)-pbinom(i,n,0.5)
if (conf>1-alpha){conff<-c(conff,i,j,conf)}
}
}
conff
length(conff)
min<-103-38;min
c<-seq(1,(length(conff)-1),3);c
for(i in c){
col<-c(rstu[conff[i]],rstu[conff[i+1]],conff[i+2])
min1<-rstu[conff[i+1]]-rstu[conff[i]]
if (min1for(k in 1:n){
conf=pbinom(n-k,n,0.5)-pbinom(k,n,0.5)
if (conf<1-alpha){loc=k-1;break}
}
print(loc)
（剩余的例题参考程序在课本）
3.6正态记分检验
例2.18
baby1<-c(4,6,9,15,31,33,36,65,77,88)
baby=(baby1-34);baby
baby.mean=mean(baby);baby.mean
例2.18
qiuzhi<-function(x){
n=length(x)
a=rep(2,n)
for (i in 1:n){
a[i]=sum(x<=x[i])
}
a
}
fuhao<-function(x,y){
n=length(x)
sgn=rep(2,n)
for(i in 1:n){
if (x[i]>y)
sgn[i]=1
else if (x[i]==y)
sgn[i]=0
else
sgn[i]=-1
}
sgn
}
n1<-length(baby)
babyzhi=qiuzhi(baby)
q=(n1+1+babyzhi)/(2*n1+2)
babysgn<-fuhao(baby,34)
babysgn=sign(baby1-34);babysgn
s=qnorm(q,0,1)
W<-t(s)%*%babysgn;W
sd<-sum((s*babysgn)^2);sd
T=W/sd;T

2.7分布的一致性检验
例2.19
shuju1<-data.frame(month=c(1:6),
customers=c(27,18,15,24,36,30));shuju1
attach(shuju1)
n<-sum(customers);n
expect<-rep(1,6)*(1/6)*n;expect

x.squ=sum((customers-expect)^2)/25;x.squ
#方法一
value<-qchisq(1-0.05,length(customers)-1);value
#方法二
pvalue<-1-pchisq(x.squ,length(customers)-1);pvalue

例2.20
shuju2<-data.frame(chongshu=c(0:6),
zhushu=c(10,24,10,4,1,0,1));shuju2
attach(shuju2)
n=sum(zhushu);n
lamda<-sum(chongshu*zhushu)/nlamda
p<-dpois(chongshu,lamda);p
n*p
x.squ=sum((zhushu^2)/(n*p))-n;x.squ
#方法一
value<-qchisq(1-0.05,length(zhushu)-1);value
#方法二
pvalue<-1-pchisq(x.squ,length(zhushu)-1);pvalue

例2.21
shuju3<-c(36,36,37,38,40,42,43,43,44,45,48,48,
50,50,51,52,53,54,54,56,57,57,57,58,58,58,58,
58,59,60,6161,61,62,62,63,63,65,66,68,68,70,
73,73,75);shuju3
n=length(shuju3)
n0=sum(shuju3<30);n0
n1=sum(shuju3>30 & shuju3<=40);n1
n2=sum(shuju3>40 & shuju3<=50);n2
n3=sum(shuju3>50 & shuju3<=60);n3
n4=sum(shuju3>60 & shuju3<=70);n4
n5=sum(shuju3>70 & shuju3<=80);n5
n6=sum(shuju3>80);n6
nn<-c(n0,n1,n2,n3,n4,n5,n6);nn #计算45位学生体重分类的频数！
shuju3.mean=mean(shuju3);shuju3.mean
shuju3.var=var(shuju3);shuju3.var
shuju3.sd=sd(shuju3);shuju3.sd
e0=pnorm(30,shuju3.mean,shuju3.sd)
e1=pnorm(40,shuju3.mean,shuju3.sd)-pnorm(30,shuju3.mean,shuju3.sd)
e2=pnorm(50,shuju3.mean,shuju3.sd)-pnorm(40,shuju3mean,shuju3.sd)
e3=pnorm(60,shuju3.mean,shuju3.sd)-pnorm(50,shuju3.mean,shuju3.sd)
e4=pnorm(70,shuju3.mean,shuju3.sd)-pnorm(60,shuju3.mean,shuju3.sd)
e5=pnorm(80,shuju3.mean,shuju3.sd)-pnorm(70,shuju3.mean,shuju3.sd)
e6=1-pnorm(80,shuju3.mean,shuju3.sd)
e=c(e0,e1,e2,e3,e4,e5,e6);e
ee=n*c(e0,e1,e2,e3,e4,e5,e6);ee
x.squ=sum((nn^2)/(ee))-n;x.squ
#方法一
value<-qchisq(1-0.05,length(ee)-1);value
#方法二
pvalue<-1-pchisq(x.squ,length(ee)-1);pvalue

例2.22
healthy<-c(87,77,92,68,80,78,84,77,81,80,80,77,92,86,
76,80,81,75,77,72,81,90,84,86,80,68,77,87,76,77,78,92,
75,80,78);healthy
ks.test(healthy,pnorm,80,6)
第三章
#Brown_Mood中位数
#Brown-Mood中位数检验程序
BM.test<-function(x,y,alt){
xy<-c(x,y)
md.xy<-median(xy) #利用中位数的检验
#md.xy<-quantile(xy,0.25) #利用p分位数的检验
t<-sum(xy>md.xy)
lx<-length(x)
ly<-length(y)
lxy<-lx+ly
A<-sum(x>md.xy)
if (alt=="greater")
{w<-1-phyper(A,lx,ly,t)}
else if (alt=="less")
{w<-phyper(A,lx,ly,t)}
conting.table=matrix(c(A,lx-A,lx,t-A,ly-(t-A),ly,t,lxy-t,lxy),3,3)
col.name<-c("X","Y","X+Y")
row.name<-c(">MXY","able,p.vlue=w)
}
例3.2
X<-c(698,688,675,656,655,648,640,639,620)
Y<-c(780,754,740,712,693,680,621)
#方法一：
BM.test(X,Y,"less")
#方法二：
XY<-c(X,Y)
md.xy<-median(XY)
t<-sum(XY>md.xy)
lx<-length(X)
ly<-length(Y)
lxy<-lx+ly
A<-sum(X>md.xy)
#没有修正时的情形
pvalue1<-pnorm(A,lx*t/(lx+ly),
sqrt(lx*ly*t*(lx+ly-t)/(lx+ly)^3));pvalue1
#修正时的情形
pvalue2<-pnorm(A,lx*t/(lx+ly)-0.5,
sqrt(lx*ly*t*(lx+ly-t)/(lx+ly)^3));pvalue2
3.2、Wilcoxon-Mann-Whitney秩和检验
#求两样本分别的秩和的程序.
Qiuzhi<-function(x,y){
n1<-length(y)
yy<-c(x,y)
wm=0
for(i in 1:n1){
wm=wm+sum(y[i]>yy,1)
}
wm
}
例3.3
weight.lowc(134,146,104,119,124,161,
107,83,113,129,97,123)
m=length(weight.low)
weight.high=c(70,118,101,85,112,132,94)
n=length(weight.high)
#方法一：
wy<-Qiuzhi(weight.low,weight.high)##wy=50
wxy<-wy-n*(n+1)/2;wxy#=22
mean<-m*n/2
var<-m*n*(m+n+1)/12
pvalue<-1-2*pnorm(wxy,mean-0.5,var);pvalue

#方法二
wilcox.test(weight.high,weight.low)

例3.4
Mx-My的R参考程序：
x1<-c(140,147,153,160,165,170,171,193)
x2<-c(130,135,138,144,148,155,168)
n1<-length(x1)
n2<-length(x2)
th.hat<-median(x2)-median(x1)
B=10000
Tboot=c(rep(0,1000)) #vector of length Bootstrap
for (i in 1:B)
{
xx1=sample(x1,5,T) #sample of size n1 with replacement from x1
xx2=sample(x2,5,T) #sample of size n2 with replacement from x2
Tboot[i]=median(xx2)-median(xx1)
}
th<-median(Tboot);th
se=sd(Tboot)
Normal.conf=c(th+qnorm(0.025,0,1)*se,th-qnorm(0.025,0,1)*se);Normal.conf
Percentile.conf=c(2*th-quantile(Tboot,0.975),2*th-quantile

(Tboot,0.025));Percentile.conf
Provotal.conf=c(quantile(Tboot,0.025),quantile(Tboot,0.975));Provotal.conf
th.hat
3.3、Mood方差检验
qiuzhi<-function(x,y){
xy<-c(x,y)
zhi<-NULL
for (i in 1:length(x)){
zhi<-c(zhi,sum(x[i]>=xy))
}
zhi
}
引例：
x1<-c(48,56,59,61,84,8791,95)
x2<-c(2,22,49,78,85,89,93,97)
zhi_x1=qiuzhi(x1,x2);zhi_x1
#zhi_x2=qiuzhi(x2,x1);zhi_x2
#var_x1=var(x1);var_x1
#var_x2=var(x2);var_x2
m=length(x1);m
n=length(x2);n
mean_R=(m+n+1)/2;mean_R
mean1=m*(m+n+1)*(m+n-1)/12;mean1
var1=m*n*(m+n+1)*(m+n+2)*(m+n-2)/180;var1
M1=sum((zhi_x1-meanR)^2);M1
p_value=2*pnorm(M1,mean1-0.5,sqrt(var1))
p_value

例3.5
X<-c(4.5,6.5,7,10,12)
Y<-c(6,7.2,8,9,9.8)
zhi_X=qiuzhi(X,Y);zhi_X
m=length(X);m
n=length(Y);n
mean_R=(m+n+1)/2;mean_R
mean2=m*(m+n+1)*(m+n-1)/12;mean2
var2=m*n*(m+n+1)*(m+n+2)*(m+n-2)/180;var2
M2=sum((zhi_X-mean_R)^2);M2
方法一：查附表9
#方法二：
p_value=2*(1-pnorm(M2,mean2-0.5,sqrt(var2)))
p_value
#方法三
Z=1/(sqrt(var2))*(M2-mean2+0.5);Z
3.4、Moses方差检验
qiuzhi<-function(x,y){
xy<-c(x,y)
zhi<-NULL
for (i in 1:length(x)){
zhi<-c(zhi,sum(x[i]>=xy))
}
zhi
}
例3.6
x1<-c(8.2,10.7,7.5,14.6,6.3,9.2,11.9,
5.6,12.8,5.2,4.9,13.5)
m1=length(x1);m1
x2<-c(4.7,6.3,5.2,6.8,5.6,4.2,
6.0,7.4,8.1,6.5)
m2=length(x2);m2
A<-matrix(x1,ncol=3);A#随机分组
a1=sample(x1,3,F)
xx2=NULL
for(i in 1:m1){
if(sum(a1==x1[i])==0) xx2=c(xx2,x1[i])
}
a2=sample(xx2,3,F)

xx3=NULL
for(i in 1:(m1-3)){
if(sum(a2==xx2[i])==0) xx3=c(xx3,x1[i])
}
a3=sample(xx3,3,F)
x11=sum((A[1,]-mean(x1))^2);x11
x12=sum((A[2,]-mean(x1))^2);x12
x13=sum((A[3,]-mean(x1))^2);x13
x14=sum((A[4,]-mean(x1))^2);x14
SSA<-c(x11,x12,x13,x14);SSA
B<-matrix(x2[1:9],ncol=3);B
y11=sum((B[1,]-mean(x2))^2);y11
y12=sum((B[2,]-mean(x2))^2);y12
y13=sum((B[3,]-mean(x2))^2);y13
SSB<-c(y11,y12,y13);SSB
zhi_SSA=qiuzhi(SSA,SSB);zhi_SSA
zhi_SSB=qiuzhi(SSB,SSA);zhi_SSB
S=sum(zhi_SSA);S
TM=S-4*(4+1)/2;TM
#方法一(查附表4)
拒绝域C=(TMW(0.975,m1,m2))
其中W(0.975,m1,m2)=m1*m2-W(0.025,m1,m2).
#方法二（Wilcoxon秩和检验）
wilcox.test(SSA,SSB)
#方法二（Mann-Whitney秩和检验）
m=length(SSA);m
n=length(SSB);n
mean_AB=m*n/2;mean_AB
var_AB=m*n*(m+n+1)/12;var_AB
p_value=1-pnorm(S,mean_AB,sqrt(var_AB));p_value
第四章
4.1、试验设计和方差分析的基本概念回顾
#R软件中单因素方差分析的函数
例4.1
#方法一：
****Analysis of Variance Model ****
y<-c(2.0,1.4,2.0,2.8,2.4,1.9,1.8,2.5,20,1.5,2.1,2.2);y
lever<-c("B","A","C","C","B","A","B","C","A","A","C","B")
x<-factor(lever);x
xy<-data.frame(y,x)
attach(xy)
aov(formula=y~x,data=xy)
aov.xy<-aov(formula=y~x,data=xy)
summary(aov.xy)

#方法二：
x1<-c(1.4,1.9,2.0,1.5)
x2<-c(2.0,2.4,1.8,2.2)
x3<-c(2.6,2.8,2.5,2.1)
y<-c(x1,x2,x3);y
y.mean<-mean(y);y.mean
ssT<-sum((y-y.mean)^2);ssT #计算总的平方和
x1.mean<-mean(x1)
x2.mean<-mean(x2)
x3.mean<-mean(x3)
sse<-sum(sum((x1-x1.mean)^2),
sum((x2-x2.mean)^2),sum((x3-x3.mean)^2));sse #计算误差平方和
sst<-ssT-sse;sst #计算组间平方和
F<-(sst/2)/(sse/(length(y)-3));F #计算方差分析的F检验统计量
#临界值的计算

value<-qf(0.95,2,length(y)-3);value
#计算p-value值
p.value<-1-pf(8,2,length(y)-3);p.value
表4.5
xueye<-c(8.4,9.4,9.8,12.2,
10.8,15.2,9.8,14.4,8.6,9.8,10.2,9.8,
8.8,9.8,8.9,12.0,8.4,9.2,8.5,9.5);xueye
sst1<-sum((xueye-mean(xueye))^2);sst1
a=matrix(xueye,ncol=5);a
quzu<-apply(a,2,sum);quzu
chuli<-apply(a,1,sum);chuli
k=5
b=4
ssb=1/4*sum(quzu^2)-sum(quzu)^2/(k*b);ssb
sst=1/5*sum(chuli^2)-sum(chuli)^2/(k*b);sst
sse=sst1-ssb-sst;sse
mssb=ssb/(k-1);mssb
msst=sst/(b-1);msst
msse=sse/(k*b-k-b+1);msse
F1=mssb/msse;F1
F2=msst/msse;F2
value1=qf(1-0.05,k-1,k*b-k-b+1)
value2=qf(1-0.05,b-1,k*b-k-b+1)

例4.3
qiuzhi<-function(w,x,y,z){
xy<-c(w,x,y,z)
zhi<-NULL
for (i in 1:length(w)){
zhi<-c(zhi,sum(w[i]>=xy))
}
zhi
}
a<-c(80,203,236,252,284,368,457,393)
b<-c(133,180,100,160)
c<-c(156,295,320,448,465,481,279)
d<-c(194,214,272,330,386,475)
azhi=qiuzhi(a,b,c,d);azhi
bzhi=qiuzhi(b,a,c,d);bzhi
czhi=qiuzhi(c,a,b,d);czhi
dzhi=qiuzhi(d,a,b,c);dzhi
H=12/(n*(n+1))*(sum(azhi)^2/length(a)+sum(bzhi)^2/length(b)
+sum(czhi)^2/length(c)+sum(dzhi)^2/length(d))-(3*(n+1))
方法一：value=qchisq(1-0.05,3);value
方法二：pvalue=1-pchisq(H,3);pvalue
mean=c(mean(a),mean(b),mean(c),mean(d))
#两两比较的程序
bjiao=function(azhi,bzhi,czhi,dzhi){
{n=length(c(azhi,bzhi,czhi,dzhi))
av=sum(azhi)/length(azhi)
bv=sum(bzhi)/length(bzhi)
se=sqrt(n*(n+1)/12*(1/length(azhi)+1/length(bzhi)))
d=abs(av-bv)
dab=d/se
huizong=c(d,se,dab,qnorm1-0.05,0,1))}
huizong
}
bjiao(azhi,bzhi,czhi,dzhi)
bjiao(czhi,dzhi,azhi,bzhi)
4.3、Jonckheere-Terpstra检验
例4.5
x=c(125,136,116,101,105,109)
y=c(122,114,131,120,119,127)
z=c(128,142,128,134,135,131,140,129)
xm=mean(x);xm
ym=mean(y);ym
zm=mean(z);zm
g=c(rep(1,6),rep(2,6),rep(3,8))
ta