完整版非参数统计第二版习题R程序Word格式.docx
《完整版非参数统计第二版习题R程序Word格式.docx》由会员分享，可在线阅读，更多相关《完整版非参数统计第二版习题R程序Word格式.docx（28页珍藏版）》请在冰豆网上搜索。
t1<
for（iin1:
length（ss））{
if（ss[i]<
0）t<
-t+1#求小于80的个数
elset1<
-t1+1求大于80的个数
}
t;
t1
t;
[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]<
rain[i+5]）t1<
-t1+1
k<
-0:
t1-1
sum（dbinom（k,5,0.5））#=0.1875
y<
-6/（2A5）;
y#=0.1875
#
（2）该地区前32年降雨量是否变化?
例2.9
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.3
3476（课本为ru=17）
t=0
16）（
rain[i+16]）t<
-t+1
}t
k1<
min（t,16-t）-1
sum（dbinom（k1,16,0.5））#=0.0002593994pbinom（max（k1）,16,0.5）#=0.0002593994y1<
-（1+16）/（2A16）;
y1#=0.0002593994plot（year,rain）
abline（v=（1971+2002）/2,col=2）lines（year,rain）
anova（lm（rain~（year）））
随机游程检验（P50）例2.8
client<
-c（"
F"
"
M"
"
clientn<
-length（client）;
nn1<
-sum（client=="
n1
n0<
-n-n1;
n0t1<
-0for（iin1:
（length（client）-1））（
if（client[i]==client[i+1]）t1<
-t1
-t1+1}
R<
-t1+1;
R#=12#findrejectionregion（不写）shuju39<
-data.frame（read.table
（"
SHUJU39.txt"
header=TRUE））;
shuju39
attach（shuju39）
sum.a=0
sum.b=0
sum.c=0
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
&
（（chanliang[i]-ma）>
0））fuhao[i]<
-"
+"
（（chanliang[i]-mb）>
C"
（（chanliang[i]-mc）>
elsefuhao[i]<
#利用上题编程解决检验的随机性
n<
-length（fuhao）;
n
n1<
-sum（fuhao=="
n0
（length（fuhao）-1））（
if（fuhao[i]==fuhao[i+1]）t1<
-t1+1;
R
#findrejectionregion
-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/（2Ai）#densityofWilcoxdistributwhen
size=N
N=19#samplesizeofexpecteddistributionofW+
-dwilxonfun（N）;
y
#计算P（W+=x）中的x取值的R参考程序！
！
p=a/（2Ai）#densityofWilcoxdistributwhensize=N
a
length（y）-1
hist（y,freq=FALSE）
lines（density（y）,col="
）
例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
（length（a）-1））{
for（jin（i+1）:
length（a））{
-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-1））{
n）{
conf=pbinom（j,n,0.5）-pbinom（i,n,0.5）
if（conf>
1-alpha）{conff<
-c（conff,i,j,conf）}
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<
min）{min<
-min1;
l<
-i}
print（col）
col1<
c（rstu[conff[l]],rstu[conff[l+1]],conff[l+2]）;
col1
例2.14“
conf=pbinom（n,n,0.5）-pbinom（0,n,0.5）;
conf
for（kin1:
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
qiuzhi<
-function（x）{
n=length（x）
a=rep（2,n）
a[i]=sum（x<
=x[i]）
-function（x,y）{
sgn=rep（2,n）
n）（
if（x[i]>
y）
sgn[i]=1
elseif（x[i]==y）
sgn[i]=0
else
sgn[i]=-1
sgn
-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）A2）;
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）
-sum（customers）;
expect<
-rep（1,6）*（1/6）*n;
expect
x.squ=sum（（customers-expect）A2）/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:
zhushu=c（10,24,10,4,1,0,1））;
shuju2
attach（shuju2）
n=sum（zhushu）;
lamda<
-sum（chongshu*zhushu）/n;
lamda
-dpois（chongshu,lamda）;
n*p
x.squ=sum（（zhushuA2）/（n*p））-n;
-qchisq（1-0.05,length（zhushu）-1）;
-1-pchisq（x.squ,length（zhushu）-1）;
例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）;
n1=sum（shuju3>
30&
shuju3<
=40）;
n2=sum（shuju3>
40&
=50）;
n2
n3=sum（shuju3>
50&
=60）;
n3
n4=sum（shuju3>
60&
=70）;
n4
n5=sum（shuju3>
70&
=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（（nnA2）/（ee））-n;
-qchisq（1-0.05,length（ee）-1）;
-1-pchisq（x.squ,length（ee）-1）;
例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,7
8,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<
-quantile（xy,0.25）#利用p分位数的检
验
-sum（xy>
md.xy）
lx<
-length（x）
ly<
-length（y）
lxy<
-lx+ly
A<
-sum（x>
if（alt=="
greater"
（w<
-1-phyper（A,lx,ly,t）}
elseif（alt=="
less"
-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<
X”,"
Y”,"
X+Y"
row.name<
MXY"
<
TOTAL"
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）
-sum（XY>
-length（X）
-length（Y）
-sum（X>
#没有修正时的情形
pvalue1<
-pnorm（A,lx*t/（lx+ly）,
-median（xy）#利用中位数的检验
Mx-My的R参考程序:
sqrt（lx*ly*t*（lx+ly-t）/（lx+ly）A3））;
pvalue1
#修正时的情形
pvalue2<
-pnorm（A,lx*t/（lx+ly）-0.5,
pvalue2
3.2、Wilcoxon-Mann-Whitney秩和检验
#求两样本分别的秩和的程序.
Qiuzhi<
-function（x,y）（
yy<
wm=0
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
-1-2*pnorm（wxy,mean-0.5,var）;
wilcox.test（weight.high,weight.low）
例3.4
x1<
-c（140,147,153,160,165,170,171,193）
x2<
-c（130,135,138,144,148,155,168）
-length（x1）
n2<
-length（x2）
th.hat<
-median（x2）-median（x1）
B=10000
Tboot=c（rep（0,1000））#vectoroflengthBootstrap
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
升级会员 