非参数统计第二版习题R程序.docx
- 文档编号:29990652
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:27
- 大小:92.78KB
非参数统计第二版习题R程序.docx
《非参数统计第二版习题R程序.docx》由会员分享,可在线阅读,更多相关《非参数统计第二版习题R程序.docx(27页珍藏版)》请在冰豆网上搜索。
非参数统计第二版习题R程序
非参数统计(第二版)习题R程序
LT
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(iin1:
length(ss)){
if(ss[i]<0)t<-t+1#求小于80的个数
elset1<-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
((chanliang[i]-mb)>0))fuhao[i]<-"+"
elseif(pinzhong[i]=="C"&((chanliang[i]-mc)>0))fuhao[i]<-"+"
elsefuhao[i]<-"-"
}
fuhao
#利用上题编程解决检验的随机性
n<-length(fuhao);n
n1<-sum(fuhao=="+");n1
n0<-n-n1;n0
t1<-0
for(iin1:
(length(fuhao)-1)){
if(fuhao[i]==fuhao[i+1])t1<-t1
elset1<-t1+1
}
R<-t1+1;R
#findrejectionregion
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
例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)#whenn=1frequencyofW+=1oro
n=1
pp=NULL#distributeofallsizefrom2toN
aa=NULL#frequencyofallsizefrom2toN
for(iin2:
N){
t=c(rep(0,i),a)
a=c(a,rep(0,i))+t
p=a/(2^i)#densityofWilcoxdistributwhensize=N
}
p
}
N=19#samplesizeofexpecteddistributionofW+
y<-dwilxonfun(N);y
#计算P(W+=x)中的x取值的R参考程序!
!
dwilxonfun=function(N){
a=c(1,1)#whenn=1frequencyofW+=1oro
n=1
pp=NULL#distributeofallsizefrom2toN
aa=NULL#frequencyofallsizefrom2toN
for(iin2:
N){
t=c(rep(0,i),a)
a=c(a,rep(0,i))+t
p=a/(2^i)#densityofWilcoxdistributwhensize=N
}
a
}
N=19#samplesizeofexpecteddistributionofW+
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(iin1:
(length(a)-1)){
for(jin(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(iin1:
(n-1)){
for(jin(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(iinc){
col<-c(rstu[conff[i]],rstu[conff[i+1]],conff[i+2])
min1<-rstu[conff[i+1]]-rstu[conff[i]]
if(min1 print(col) } col1<-c(rstu[conff[l]],rstu[conff[l+1]],conff[l+2]);col1 min 例2.14“ stu<-c(82,53,70,73,103,71,69, 80,54,38,87,91,62,75,65,77);stu alpha=0.05 n=length(stu);n conf=pbinom(n,n,0.5)-pbinom(0,n,0.5);conf for(kin1: 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(iin1: n){ a[i]=sum(x<=x[i]) } a } fuhao<-function(x,y){ n=length(x) sgn=rep(2,n) for(iin1: n){ if(x[i]>y) sgn[i]=1 elseif(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)/n;lamda 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,61,61,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,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=="greater") {w<-1-phyper(A,lx,ly,t)} elseif(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"," 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,"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(iin1: 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))#vectoroflengthBootstrap for(iin1: B) { xx1=sample(x1,5,T)#sampleofsizen1withreplacementfromx1 xx2=sample(x2,5,T)#sampleofsizen2withreplacementfromx2 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(iin1: 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(iin1: 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<-mat
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 参数 统计 第二 习题 程序