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,shuju3.mean,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==\ {w<-1-phyper(A,lx,ly,t)} else if (alt==\ {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(\
row.name<-c(\
dimnames(conting.table)<-list(row.name, col.name) list(contingency.table=conting.table,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,\#方法二： 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.low=c(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,87,91,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-mean_R)^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=(TM
其中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,2.0,1.5,2.1,2.2);y lever<-c(\
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,qnorm(1-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)) tapply(c(x,y,z),g,median) JT.test(data=t(c(x,y,z)),class=g) Wij<-function(x,y){