统计建模与R软件课后答案.docx

统计建模与R软件课后答案.docx

《统计建模与R软件课后答案.docx》由会员分享，可在线阅读，更多相关《统计建模与R软件课后答案.docx（32页珍藏版）》请在冰豆网上搜索。

统计建模与R软件课后答案

第二章

2.1

>x<-c（1,2,3）;y<-c（4,5,6）

>e<-c（1,1,1）

>z<-2*x+y+e;z

[1]71013

>z1<-crossprod（x,y）;z1

[,1]

[1,]32

>z2<-outer（x,y）;z2

[,1][,2][,3]

[1,]456

[2,]81012

[3,]121518

2.2

（1）>A<-matrix（1:

20,nrow=4）;B<-matrix（1:

20,nrow=4,byrow=T）

>C<-A+B;C

（2）>D<-A%*%B;D

（3）>E<-A*B;E

（4）>F<-A[1:

3,1:

3]

（5）>G<-B[,-3]

2.3

>x<-c（rep（1,5）,rep（2,3）,rep（3,4）,rep（4,2））;x

2.4

>H<-matrix（nrow=5,ncol=5）

>for（iin1:

5）

+for（jin1:

5）

+H[i,j]<-1/（i+j-1）

（1）>det（H）

（2）>solve（H）

（3）>eigen（H）

2.5

>studentdata<-data.frame（姓名=c（'张三','李四','王五','赵六','丁一'）

+,性别=c（'女','男','女','男','女'）,年龄=c（'14','15','16','14','15'）,

+身高=c（'156','165','157','162','159'）,体重=c（'42','49','41.5','52','45.5'））

2.6

>write.table（studentdata,file='student.txt'）

>write.csv（studentdata,file='student.csv'）

2.7

count<-function（n）

{

if（n<=0）

print（'要求输入一个正整数'）

else{

repeat{

if（n%%2==0）

n<-n/2

else

n<-（3*n+1）

if（n==1）break

}

print（'运算成功'）}

}

第三章

3.1

首先将数据录入为x。

利用data_outline函数。

如下

>data_outline（x）

3.2

>hist（x,freq=F）

>lines（density（x）,col='red'）

>y<-min（x）:

max（x）

>lines（y,dnorm（y,73.668,3.9389）,col='blue'）

>plot（ecdf（x）,verticals=T,do.p=F）

>lines（y,pnorm（y,73.668,3.9389））

>qqnorm（x）

>qqline（x）

3.3

>stem（x）

>boxplot（x）

>fivenum（x）

3.4

>shapiro.test（x）

>ks.test（x,'pnorm',73.668,3.9389）

One-sampleKolmogorov-Smirnovtest

data:

x

D=0.073,p-value=0.6611

alternativehypothesis:

two-sided

Warningmessage:

Inks.test（x,"pnorm",73.668,3.9389）:

tiesshouldnotbepresentfortheKolmogorov-Smirnovtest

这里出现警告信息是因为ks检验要求样本数据是连续的，不允许出现重复值

3.5

>x1<-c（2,4,3,2,4,7,7,2,2,5,4）;x2<-c（5,6,8,5,10,7,12,12,6,6）;x3<-c（7,11,6,6,7,9,5,5,10,6,3,10）

>boxplot（x1,x2,x3,names=c（'x1','x2','x3'）,vcol=c（2,3,4））

>windows（）

>plot（factor（c（rep（1,length（x1））,rep（2,length（x2））,rep（3,length（x3））））,c（x1,x2,x3））

3.6

>rubber<-data.frame（x1=c（65,70,70,69,66,67,68,72,66,68）,

+x2=c（45,45,48,46,50,46,47,43,47,48）,x3=c（27.6,30.7,31.8,32.6,31.0,31.3,37.0,33.6,33.1,34.2））

>plot（rubber）

具体有相关关系的两个变量的散点图要么是从左下角到右上角（正相关），要么是从左上角到右下角（负相关）。

从上图可知所有的图中偶读没有这样的趋势，故均不相关。

3.7

（1）>student<-read.csv（'3.7.csv'）

>attach（student）

>plot（体重~身高）

（2）>coplot（体重~身高|性别）

（3）>coplot（体重~身高|年龄）

（4）>coplot（体重~身高|年龄+性别）

只列出（4）的结果，如下图

3.8

>x<-seq（-2,3,0.5）;y<-seq（-1,7,0.5）

>f<-function（x,y）

+x^4-2*x^2*y+x^2-2*x*y+2*y^2+9*x/2-4*y+4

>z<-outer（x,y,f）

>contour（x,y,z,levels=c（0,1,2,3,4,5,10,15,20,30,40,50,60,80,100）,col='blue'）

>windows（）

>persp（x,y,z,theta=30,phi=30,expand=0.7,col='red'）

3.9

>cor.test（身高,体重）

根据得出的结果看是相关的。

具体结果不再列出

3.10

>df<-read.csv（'48名求职者得分.csv'）

>stars（df）

然后按照G的标准来画出星图

>attach（df）

>df$G1<-（SC+LC+SMS+DRV+AMB+GSP+POT）/7

>df$G2<-（FL+EXP+SUIT）/3

>df$G3<-（LA+HON+KJ）/3

>df$G4<-AA

>df$G5<-APP

>a<-scale（df[,17:

21]）

>stars（a）

这里从17开始取，是因为在df中将ID也作为了一列

3.11

使用P159已经编好的函数unison，接着上题，直接有

>unison（a）

第四章

4.1

（1）先求矩估计。

总体的期望为

。

因此我们有

。

可解得a=（2*E（

）-1）/（1-E（

））.因此我们用样本的均值来估计a即可。

在R中实现如下

>x<-c（0.1,0.2,0.9,0.8,0.7,0.7）

>（2*mean（x）-1）/（1-mean（x））

[1]0.3076923

（2）采用极大似然估计

首先求出极大似然函数为

再取对数为

最后求导

好了下面开始用R编程求解，注意此题中n=6.

方法一、

使用unniroot函数

>f<-function（a）6/（a+1）+sum（log（x））

>uniroot（f,c（0,1））

方法二、

使用optimize函数

>g<-function（a）6*log（a+1）+a*sum（log（x））

>optimize（g,c（0,1）,maximum=T）

4.2

用极大似然估计得出

.现用R求解如下

>x<-c（rep（5,365）,rep（15,245）,rep（25,150）,rep（35,100）,rep（45,70）,rep（55,45）,rep（65,25））

>1000/sum（x）

4.3

换句话讲，就是用该样本来估计泊松分布中的参数，然后求出该分布的均值。

我们知道泊松分布中的参数

，既是均值又是方差。

因此我们只需要用样本均值作矩估计即可

在R中实现如下

>x<-c（rep（0,17）,rep（1,20）,rep（2,10）,rep（3,2）,rep（4,1））

>mean（x）

[1]1

4.4

>f<-function（x）{

+obj<-c（-13+x[1]+（（5-x[2]）*x[2]-2）*x[2],（-29+x[1]+（（x[2]+1）*x[2]-14）*x[2]））

+sum（obj^2）}

>nlm（f,c（0.5,-2））

4.5

在矩估计中，正态分布总体的均值用样本的均值估计。

故在R中实现如下

>x<-c（54,67,68,78,70,66,67,70,65,69）

>mean（x）

[1]67.4

然后用t.test作区间估计，如下

>t.test（x）

>t.test（x,alternative='less'）

>t.test（x,alternative='greater'）

此时我们只需要区间估计的结果，所以我们只看t.test中的关于置信区间的输出即可。

t.test同时也给出均值检验的结果，但是默认mu=0

并不是我们想要的。

下面我们来做是否低于72的均值假设检验。

如下

>t.test（x,alternative='greater',mu=72）

OneSamplet-test

data:

x

t=-2.4534,df=9,p-value=0.9817

alternativehypothesis:

truemeanisgreaterthan72

95percentconfidenceinterval:

63.96295Inf

sampleestimates:

meanofx

67.4

结果说明：

我们的备择假设是比72要大，但是p值为0.9817，所以我们不接受备择假设，接受原假设比72小。

因此这10名患者的平均脉搏次数比正常人要小。

4.6

我们可以用两种方式来做一做

>x<-c（140,137,136,140,145,148,140,135,144,141）

>y<-c（135,118,115,140,128,131,130,115,131,125）

>t.test（x,y,var.equal=T）

>t.test（x-y）

结果不再列出，但是可以发现用均值差估计和配对数据估计的结果的数值有一点小小的差别。

但得出的结论是不影响的（他们的期望差别很大）

4.7

>A<-c（0.143,0.142,0.143,0.137）

>B<-c（0.140,0.142,0.136,0.138,0.140）

>t.test（A,B）

4.8

>x<-c（140,137,136,140,145,148,140,135,144,141）

>y<-c（135,118,115,140,128,131,130,115,131,125）

>var.test（x,y）

>t.test（x,y,var.equal=F）

4.9

泊松分布的参数就等于它的均值也等于方差。

我们直接用样本均值来估计参数即可，然后作样本均值0.95的置信区间即可。

>x<-c（rep（0,7）,rep（1,10）,rep（2,12）,rep（3,8）,rep（4,3）,rep（5,2））

>mean（x）

[1]1.904762

>t.test（x）

4.10

正态总体均值用样本均值来估计。

故如下

>x<-c（1067,919,1196,785,1126,936,918,1156,920,948）

>t.test（x,alternative='greater'）

注意greater才是求区间下限的（都比它大的意思嘛）

第五章

5.1

这是一个假设检验问题，即检验油漆作业工人的血小板的均值是否为225.在R中实现如下

>x<-scan（）

1:

220188162230145160238188247113

11:

126245164231256183190158224175

21:

Read20items

>t.test（x,mu=225）

5.2

考察正态密度函数的概率在R中的计算。

首先我们要把该正态分布的均值和方差给估计出来，这个就利用样本即可。

然后用pnorm函数来计算大于1000的概率。

如下

>x<-c（1067,919,1196,785,1126,936,918,1156,920,948）

>pnorm（1000,mean（x）,sd（x））

[1]0.5087941

>1-0.5087941

[1]0.4912059

5.3

这是检验两个总体是否存在差异的问题。

可用符号检验和wilcoxon秩检验。

两种方法实现如下

>x<-c（113,120,138,120,100,118,138,123）

>y<-c（138,116,125,136,110,132,130,110）

>binom.test（sum（x

p-value=1

>wilcox.test（x,y,exact=F）

p-value=0.792

可见无论哪种方法P值都大于0.05，故接受原假设，他们无差异

5.4

（1）采用w检验法

>x<-c（-0.7,-5.6,2,2.8,0.7,3.5,4,5.8,7.1,-0.5,2.5,-1.6,1.7,3,0.4,4.5,4.6,2.5,6,-1.4）

>y<-c（3.7,6.5,5,5.2,0.8,0.2,0.6,3.4,6.6,-1.1,6,3.8,2,1.6,2,2.2,1.2,3.1,1.7,-2）

>shapiro.test（x）

>shapiro.test（y）

采用ks检验法

>ks.test（x,'pnorm',mean（x）,sd（x））

>ks.test（y,'pnorm',mean（y）,sd（y））

采用pearson拟合优度法对x进行检验

>A<-table（cut（x,br=c（-2,0,2,4,6,8）））

>A

（-2,0]（0,2]（2,4]（4,6]（6,8]

44641

发现A中有频数小于5，故应该重新调整分组

>A<-table（cut（x,br=c（-2,2,4,8）））

>A

（-2,2]（2,4]（4,8]

865

然后再计算理论分布

>p<-pnorm（c（-2,2,4,8）,mean（x）,sd（x））

>p<-c（p[2],p[3]-p[2],1-p[3]）

最后检验

>chisq.test（A,p=p）

采用pearson拟合优度法对y进行检验

>B<-table（cut（y,br=c（-2.1,1,2,4,7）））

>B

（-2.1,1]（1,2]（2,4]（4,7]

5555

>p<-pnorm（c（1,2,4）,mean（y）,sd（y））

>p<-c（p[1],p[2]-p[1],p[3]-p[2],1-p[3]）

>chisq.test（B,p=p）

以上的所有结果都不再列出，结论是试验组和对照组都是来自正态分布。

（2）>t.test（x,y,var.equal=F）

>t.test（x,y,var.equal=T）

>t.test（x,y,paired=T）

结论是均值无差异

（3）>var.test（x,y）

结论是方差相同

由以上结果可以看出这两种药的效果并无二致

5.5

（1）对新药组应用chisq.test检验（也可用ke.test检验）

>x<-c（126,125,136,128,123,138,142,116,110,108,115,140）

>y<-c（162,172,177,170,175,152,157,159,160,162）

>p<-pnorm（c（105,125,145）,mean（x）,sd（x））

>p<-c（p[2],1-p[2]）

>chisq.test（A,p=p）

对对照组用ks.test检验

>ks.test（y,'pnorm',mean（y）,sd（y））

结论是他们都服从正态分布

（2）>var.test（x,y）

结论是方差相同

（3）>wilcox.test（x,y,exact=F）

结果是有差别

5.6

明显是要检验二项分布的p值是否为0.147.R实现如下

>binom.test（57,400,p=0.147）

结果是支持

5.7

也就是检验二项分布中的p值是否大于0.5

>binom.test（178,328,p=0.5,alternative='greater'）

结果是不能认为能增加比例

5.8

就是检验你的样本是否符合那个分布

>chisq.test（c（315,101,108,32）,p=c（9,3,3,1）/16）

结果显示符合自由组合规律

5.9

又是检验一个总体是否符合假定分布。

>x<-0:

5;y<-c（92,68,28,11,1,0）

>z<-rep（x,y）

>A<-table（cut（z,br=c（-1,0,1,2,5）））

>q<-ppois（c（0,1,2,5）,mean（z））

>p<-c（q[1],q[2]-q[1],q[3]-q[2],1-q[3]）

>chisq.test（A,p=p）

结论是符合泊松分布

5.10

>x<-c（2.36,3.14,7.52,3.48,2.76,5.43,6.54,7.41）

>y<-c（4.38,4.25,6.53,3.28,7.21,6.55）

>ks.test（x,y）

5.11

即列联表的的独立性检验

>x<-c（358,229,2492,2754）

>dim（x）<-c（2,2）

>chisq.test（x）或>fisher.test（x）

结论是有影响

5.12

>x<-c（45,12,10,46,20,28,28,23,30,11,12,35）

>dim（x）<-c（4,3）

>chisq.test（x）

结果是相关

5.13

>x<-c（3,4,6,4）

>dim（x）<-c（2,2）

>fisher.test（x）

结果显示工艺对产品质量无影响

5.14

即检验两种研究方法是否有差异

>x<-c（58,2,3,1,42,7,8,9,17）

>dim（x）<-c（3,3）

>mcnemar.test（x,correct=F）

结果表明两种检测方法有差异

5.15

>x<-c（13.32,13.06,14.02,11.86,13.58,13.77,13.51,14.42,14.44,15.43）

>binom.test（sum（x>14.6）,length（x）,al='l'）

>wilcox.test（x,mu=14.6,al='l',exact=F）

结果表明是在中位数之下

5.16

（1）

（2）（3）

>x<-scan（）

1:

48.033.037.548.042.540.042.036.011.322.0

11:

36.027.314.232.152.038.017.320.021.046.1

21:

Read20items

>y<-scan（）

1:

37.041.023.417.031.540.031.036.05.711.5

11:

21.06.126.521.344.528.022.620.011.022.3

21:

Read20items

>binom.test（sum（x

>wilcox.test（x,y,paired=T,exact=F）

>wilcox.test（x,y,exact=F）

（4）>ks.test（x,'pnorm',mean（x）,sd（x））

>ks.test（y,'pnorm',mean（y）,sd（y））

>var.test（x,y）

由以上检验可知数据符合正态分布且方差相同，故可做t检验

>t.test（x,y）

可以发现他们的均值是有差别的

（5）综上所述，Wilcoxon符号秩检验的差异检出能力最强，符号检验的差异检出最弱。

5.17

>x<-c（24,17,20,41,52,23,46,18,15,29）

>y<-c（8,1,4,7,9,5,10,3,2,6）

>cor.test（x,y,method='spearman'）

>cor.test（x,y,method='kendall'）

有关系的

5.18

>x<-1:

5

>y<-c（rep（x,c（0,1,9,7,3）））

>z<-c（rep（x,c（2,2,11,4,1）））

>wilcox.test（y,z,exact=F）

结果显示这两种疗法没什么区别

第六章

6.1

（1）>snow<-data.frame（X=c（5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4）,

+Y=c（1907,1287,2700,2373,3260,3000,1947,2273,3113,2493））

>plot（snow$X,snow$Y）

结论是有线性关系的。

（2）（3）

>lm.sol<-lm（Y~1+X,data=snow）;summary（lm.sol）

结果是方程是显著的

（4）>predict（lm.sol,data.frame（X=7）,interval='prediction',level=0.95）

fitlwrupr

12690.2272454.9712925.484

6.2

（1）

（2）

>soil<-data.frame（X1=c（0.4,0.4,3.1,0.6,4.7,1.7,9.4,10.1,11.6,12.6,

+10.9,23.1,23.1,21.6,23.1,1.9,26.8,29.9）,X2=c（52,23,19,34,24,65,44,31,

+29,58,37,46,50,44,56,36,58,51）,X3=c（158,163,37,157,59,123,46,117,

+173,112,111,114,134,73,168,143,202,124）,Y=c（64,60,71,61,54,77,81,

+93,93,51,76,96,77,93,95,54,168,99））

>lm.sol<-lm（Y~1+X1+X2+X3,data=soil）;summary（lm.sol）

我们发现X2和X3的系数没有通过t检验。

但是整个方程通过了检验。

（3）>lm.ste