完整版非参数统计第二版习题R程序.docx
《完整版非参数统计第二版习题R程序.docx》由会员分享，可在线阅读，更多相关《完整版非参数统计第二版习题R程序.docx（28页珍藏版）》请在冰点文库上搜索。
完整版非参数统计第二版习题R程序
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
for（iin1:
5）{
if（rain[i]}
t1
k<-0:
t1-1
sum（dbinom（k,5,0.5））#=0.1875
y<-6/（2^5）;y#=0.1875
#
（2）该地区前32年降雨量是否变化？
t=0
for（iin1:
16）{
if（rain[i]}
t
k1<-0:
min（t,16-t）-1
sum（dbinom（k1,16,0.5））#=0.0002593994
pbinom（max（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（iin1:
（length（client）-1））{
if（client[i]==client[i+1]）t1<-t1
elset1<-t1+1
}
R<-t1+1;R#=12
#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#=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（iin1:
length（id））{
if（pinzhong[i]=="A"）sum.a<-sum.a+chanliang[i]
elseif（pinzhong[i]=="B"）sum.b<-sum.b+chanliang[i]
elsefuhao<-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（iin1:
length（id））{
if（pinzhong[i]=="A"&（（chanliang[i]-ma）>0））fuhao[i]<-"+"
elseif（pinzhong[i]=="B"&（（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（min1print（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-qnor
升级会员 