非全参数统计第二版习题R程序.docx

1、非全参数统计第二版习题R程序P37.例2.1build.price-c(36,32,31,25,28,36,40,32,41,26,35,35,32,87,33,35);build.pricehist(build.price,freq=FALSE)#直方图lines(density(build.price),col=red)#连线#方法一：m-mean(build.price);m#均值D-var(build.price)#方差SD-sd(build.price)#标准差St=(m-37)/(SD/sqrt(length(build.price);t#t统计量计算检验统计量t=1 -0.141

2、2332#方法二：t.test(build.price-37)#课本第38页例2.2binom.test(sum(build.price37),length(build.price),0.5)#课本40页例2.3P-2*(1-pnorm(1.96,0,1);P1 0.04999579P1 例2.4 p p1-2*(pnorm(-0.9487,0,1);p11 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(

3、scores)#输入向量求长度ss-c(scores-80);sst-0t1-0for(i in 1:length(ss)if (ssi0) t-t+1#求小于80的个数else t1 t;t11 131 15binom.test(sum(scores80),length(scores),0.75)p-value = 0.0014360.01Cox-Staut趋势存在性检验P47例2.6year-1971:2002;yearlength(year)rain-c(206,223,235,264,229,217,188,204,182,230,223,227,242,238,207,208,216

4、,233,233,274,234,227,221,214,226,228,235,237,243,240,231,210)length(rain)#(1)该地区前10年降雨量是否变化？t1=0for (i in 1:5)if (rainiraini+5) t1-t1+1t1k-0:t1-1sum(dbinom(k,5,0.5)# =0.1875y-6/(25);y# =0.1875#(2)该地区前32年降雨量是否变化？t=0for (i in 1:16)if (rainiraini+16) t-t+1tk1-0:min(t,16-t)-1sum(dbinom(k1,16,0.5)# =0.00

5、02593994pbinom(max(k1),16,0.5)#= 0.0002593994y1-(1+16)/(216);y1#=0.0002593994plot(year,rain)abline(v=(1971+2002)/2,col=2)lines(year,rain)anova(lm(rain(year)随机游程检验（P50）例2.8client-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);clientn-length(client);nn1-sum(client=M);n1n0-n-n1;n0t1-0for (i in

6、1:(length(client)-1)if (clienti=clienti+1) t1-t1else t1-t1+1R-t1+1;R#=12#find rejection region（不写）rl-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0);rlru-2*n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0);ru#=15.33476（课本为ru=17）例2.9shuju39-data.frame(read.table(SHUJU39.txt,header=TRUE);shuju39attach(shuju39)sum.a=0sum.b=0sum

7、.c=0for (i in 1:length(id)if (pinzhongi=A) sum.a-sum.a+chanliangielse if (pinzhongi=B) sum.b-sum.b+chanliangielse fuhao-sum.c-sum.c+chanliangisum.a;sum.b;sum.cma-sum.a/4mb-sum.b/4mc-sum.c/4ma;mb;mcfuhao0) fuhaoi0) fuhaoi0) fuhaoi-+else fuhaoi-fuhao#利用上题编程解决检验的随机性n-length(fuhao);nn1-sum(fuhao=+);n1n0

8、-n-n1;n0t1-0for (i in 1:(length(fuhao)-1)if (fuhaoi=fuhaoi+1) t1-t1else t1-t1+1R-t1+1;R#find rejection regionrl-1+2*n1*n0/(n1+n0)*(1-1.96/sqrt(n1+n0);rlru-2*n1*n0/(n1+n0)*(1+1.96/sqrt(n1+n0);ru例2.10（P52）library(quadprog)# 不存在叫quadprog这个名字的程辑包library(zoo)# 不存在叫zoo这个名字的程辑包library(tseries)# 不存在叫tseries

9、这个名字的程辑包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)；run1y=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 on=1pp=NULL #distr

10、ibute of all size from 2 to Naa=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/(2i) #density of Wilcox distribut when size=N pN=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 f

11、requency of W+=1 or on=1pp=NULL #distribute of all size from 2 to Naa=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/(2i) #density of Wilcox distribut when size=N aN=19 #sample size of expected distribution of W+y-dwilxonfun(N);length(y)-1hist(y,freq=FALS

12、E)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.num320);ceo.numn=length(ceo)binom.test(ceo.num,n,0.5)例2.13(P61)a-c(62,70,74,75,77,80,83,85,88)walsh-NULLfor (i in 1:(length(a)-1) for (j

13、in (i+1):length(a) walsh-c(walsh,(ai+aj)/2) walsh=c(walsh,a)NW=length(walsh);NWmedian(walsh)2.5单组数据的位置参数置信区间估计(P61)例2.14stu-c(82,53,70,73,103,71,69,80,54,38,87,91,62,75,65,77);stualpha=0.05rstu-sort(stu);rstuconff1-alpha)conff-c(conff,i,j,conf) confflength(conff)min-103-38;minc-seq(1,(length(conff)-

14、1),3);cfor(i in c)col-c(rstuconffi,rstuconffi+1,conffi+2)min1-rstuconffi+1-rstuconffiif (min1min)min-min1;l-iprint(col)col1-c(rstuconffl,rstuconffl+1,conffl+2);col1min例2.14“stu-c(82,53,70,73,103,71,69,80,54,38,87,91,62,75,65,77);stualpha=0.05n=length(stu);nconf=pbinom(n,n,0.5)-pbinom(0,n,0.5);conffo

15、r(k in 1:n) conf=pbinom(n-k,n,0.5)-pbinom(k,n,0.5) if (conf1-alpha)loc=k-1;break print(loc)（剩余的例题参考程序在课本）3.6正态记分检验例2.18baby1-c(4,6,9,15,31,33,36,65,77,88)baby=(baby1-34);babybaby.mean=mean(baby);baby.mean例2.18qiuzhi-function(x) n=length(x) a=rep(2,n) for (i in 1:n) ai=sum(x=xi) afuhaoy) sgni=1 else

16、if (xi=y) sgni=0 else sgni=-1 sgnn1-length(baby)babyzhi=qiuzhi(baby)q=(n1+1+babyzhi)/(2*n1+2)babysgn-fuhao(baby,34)babysgn=sign(baby1-34);babysgns=qnorm(q,0,1)W-t(s)%*%babysgn;Wsd-sum(s*babysgn)2);sdT=W/sd;T2.7分布的一致性检验例2.19shuju1-data.frame(month=c(1:6),customers=c(27,18,15,24,36,30);shuju1attach(sh

17、uju1)n-sum(customers);nexpect-rep(1,6)*(1/6)*n;expectx.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.20shuju2-data.frame(chongshu=c(0:6),zhushu=c(10,24,10,4,1,0,1);shuju2attach(shuju2)n=sum(zhushu);nlamda

18、-sum(chongshu*zhushu)/n;lamdap-dpois(chongshu,lamda);pn*px.squ=sum(zhushu2)/(n*p)-n;x.squ#方法一value-qchisq(1-0.05,length(zhushu)-1);value#方法二pvalue-1-pchisq(x.squ,length(zhushu)-1);pvalue例2.21shuju3-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,61,61,61,6

19、2,62,63,63,65,66,68,68,70,73,73,75);shuju3n=length(shuju3)n0=sum(shuju330 & shuju340 & shuju350 & shuju360 & shuju370 & shuju380);n6nn-c(n0,n1,n2,n3,n4,n5,n6);nn #计算45位学生体重分类的频数！shuju3.mean=mean(shuju3);shuju3.meanshuju3.var=var(shuju3);shuju3.varshuju3.sd=sd(shuju3);shuju3.sde0=pnorm(30,shuju3.mean

20、,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,s

21、huju3.sd)-pnorm(70,shuju3.mean,shuju3.sd)e6=1-pnorm(80,shuju3.mean,shuju3.sd)e=c(e0,e1,e2,e3,e4,e5,e6);eee=n*c(e0,e1,e2,e3,e4,e5,e6);eex.squ=sum(nn2)/(ee)-n;x.squ#方法一value-qchisq(1-0.05,length(ee)-1);value#方法二pvalue-1-pchisq(x.squ,length(ee)-1);pvalue例2.22healthy-c(87,77,92,68,80,78,84,77,81,80,80,7

22、7,92,86,76,80,81,75,77,72,81,90,84,86,80,68,77,87,76,77,78,92,75,80,78);healthyks.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分位数的检验 tmd.xy) lx-length(x) ly-length(y) lxy-lx+ly Amd.xy) if (alt=

23、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.nameMXY,MXY,TOTAL) dimnames(conting.table)-list(row.name, col.name) list(contingency.table=conting.table,p.vlue=w)例3.2X-c(698,688,675,656,6

24、55,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)tmd.xy)lx-length(X)ly-length(Y)lxy-lx+lyAmd.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)

25、;pvalue23.2、Wilcoxon-Mann-Whitney秩和检验#求两样本分别的秩和的程序.Qiuzhi-function(x,y)n1-length(y)yyyy,1) wm例3.3weight.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=50wxy-wy-n*(n+1)/2;wxy

26、#=22mean-m*n/2var-m*n*(m+n+1)/12pvalue-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=10000Tboot=c(rep(0,1000) #vector of length Bootstrapfor (i in 1:B)xx1=sample(x1,5,T) #sample of size n1 with replacement from x1xx2=sample(x2,5,T) #sample of size n2 with replacement from x2Tbooti=median(xx2)-median(xx1)th-median(Tboot);thse=sd(Tboot)Normal.conf=c(th+qnorm(0.025,0,1)*se,th-qnorm(