完整版非参数统计第二版习题R程序Word格式.docx

完整版非参数统计第二版习题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