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

资源描述

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

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

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

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

Themanuscriptwasrevisedontheeveningof2021

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

第二章

>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

（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:

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

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

>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）

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

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

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

>（studentdata,file=''）

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（'运算成功'）}

}

第三章

首先将数据录入为x。

利用data_outline函数。

如下

>data_outline（x）

>hist（x,freq=F）

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

>y<-min（x）:

max（x）

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

>plot（ecdf（x）,verticals=T,=F）

>lines（y,pnorm（y,,）

>qqnorm（x）

>qqline（x）

>stem（x）

>boxplot（x）

>fivenum（x）

>（x）

>（x,'pnorm',,

One-sampleKolmogorov-Smirnovtest

data:

D=,p-value=

alternativehypothesis:

two-sided

Warningmessage:

In（x,"pnorm",,:

tiesshouldnotbepresentfortheKolmogorov-Smirnovtest

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

>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））

>rubber<（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,,,,,,,,,）

>plot（rubber）

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

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

（1）

>attach（student）

>plot（体重~身高）

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

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

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

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

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

>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=,col='red'）

>（身高,体重）

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

具体结果不再列出

>df<（'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也作为了一列

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

>unison（a）

第四章

（1）先求矩估计。

总体的期望为

。

因此我们有

。

可解得a=（2*E（

）-1）/（1-E（

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

在R中实现如下

>x<-c,,,,,

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

[1]

（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）

用极大似然估计得出

.现用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）

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

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

，既是均值又是方差。

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

在R中实现如下

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

>mean（x）

[1]1

>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,-2））

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

故在R中实现如下

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

>mean（x）

[1]

然后用作区间估计，如下

>（x）

>（x,alternative='less'）

>（x,alternative='greater'）

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

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

并不是我们想要的。

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

如下

>（x,alternative='greater',mu=72）

OneSamplet-test

data:

t=,df=9,p-value=

alternativehypothesis:

truemeanisgreaterthan72

95percentconfidenceinterval:

Inf

sampleestimates:

meanofx

结果说明：

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

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

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

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

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

>（x,y,=T）

>（x-y）

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

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

>A<-c,,,

>B<-c,,,,

>（A,B）

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

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

>（x,y）

>（x,y,=F）

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

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

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

>mean（x）

[1]

>（x）

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

故如下

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

>（x,alternative='greater'）

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

第五章

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

>x<-scan（）

220188162230145160238188247113

11:

126245164231256183190158224175

21:

Read20items

>（x,mu=225）

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

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

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

如下

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

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

[1]

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

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

两种方法实现如下

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

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

>（sum（x

p-value=1

>（x,y,exact=F）

p-value=

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

（1）采用w检验法

>x<-c,,2,,,,4,,,,,,,3,,,,,6,

>y<-c,,5,,,,,,,,6,,2,,2,,,,,-2）

>（x）

>（y）

采用ks检验法

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

>（y,'pnorm',mean（y）,sd（y））

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

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

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

44641

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

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

（-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]）

最后检验

>（A,p=p）

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

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

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]）

>（B,p=p）

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

（2）>（x,y,=F）

>（x,y,=T）

>（x,y,paired=T）

结论是均值无差异

（3）>（x,y）

结论是方差相同

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

（1）对新药组应用检验（也可用检验）

>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]）

>（A,p=p）

对对照组用检验

>（y,'pnorm',mean（y）,sd（y））

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

（2）>（x,y）

结论是方差相同

（3）>（x,y,exact=F）

结果是有差别

>（57,400,p=

结果是支持

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

>（178,328,p=,alternative='greater'）

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

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

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

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

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

>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]）

>（A,p=p）

结论是符合泊松分布

>x<-c,,,,,,,

>y<-c,,,,,

>（x,y）

即列联表的的独立性检验

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

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

>（x）或>（x）

结论是有影响

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

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

>（x）

结果是相关

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

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

>（x）

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

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

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

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

>（x,correct=F）

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

>x<-c,,,,,,,,,

>（sum（x>,length（x）,al='l'）

>（x,mu=,al='l',exact=F）

结果表明是在中位数之下

（1）

（2）（3）

>x<-scan（）

11:

21:

Read20items

>y<-scan（）

11:

21:

Read20items

>（sum（x

>（x,y,paired=T,exact=F）

>（x,y,exact=F）

（4）>（x,'pnorm',mean（x）,sd（x））

>（y,'pnorm',mean（y）,sd（y））

>（x,y）

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

>（x,y）

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

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

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

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

>（x,y,method='spearman'）

>（x,y,method='kendall'）

有关系的

>x<-1:

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

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

>（y,z,exact=F）

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

第六章

（1）>snow<（X=c,,,,,,,,,,

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

>plot（snow$X,snow$Y）

结论是有线性关系的。

（2）（3）

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

结果是方程是显着的

（4）>predict,（X=7）,interval='prediction',level=

fitlwrupr

（1）

（2）

>soil<（X1=c,,,,,,,,,,

+,,,,,,,,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（Y~1+X1+X2+X3,data=soil）;summary

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

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

（3）><-step

>summary

可以发现新模型只含有X1和X3，但是X3的系数还是不显着。

接下来考虑用drop1函数处理

>drop1

发现去掉X3残差升高最小，AIC只是有少量增加。

因此应该去掉X3

><-lm（Y~X1,data=soil）;summary

此时发现新模型系数显着且方程显着

（1）>da<（X=c（1,1,1,1,2,2,2,3,3,3,4,4,4,5,6,6,6,7,7,7,8,8,8,

+9,11,12,12,12）,Y=c,,,,,,,,,,,,

+,,,,,,,,,,,,,,

+,）

>plot（da$X,da$Y）

><-lm（Y~X,data=da）

>abline

（2）>summary

全部通过

（3）>plot,1）

>windows（）

>plot,3）

可以观察到误差符合等方差的。

但是有残差异常值点24,27,28.

（4）><-update,sqrt（.）~.）

>summary

都通过检验

>plot（da$X,da$Y）

>abline

>windows（）

>plot,1）

>windows（）

>plot,3）

可以发现还是有残差离群值24,28

><-lm（Y~1+X1+X2,data=toothpaste）;summary

>plot,3）

通过函数发现5,8,9,24对样本影响较大，可能是异常值点，而通过残差图发现5是残差离群点，但是整个残差还是在[-2,2]之内的。

因此可考虑剔除5,8,9,24点再做拟合。

><-lm（Y~1+X1+X2,data=toothpaste,subset=c（-5,-8,-9,-24））

>windows（）

>plot,3）

>summary

我们发现模型的残差都控制在[,]之内，而且方程系数和方程本身也都通过检验。

>cement<（X1=c（7,1,11,11,7,11,3,1,2,21,1,11,10）,

+X2=c（26,29,56,31,52,55,71,31,54,47,40,66,68）,

+X3=c（6,15,8,8,6,9,17,22,18,4,23,9,8）,

+X4=c（60,52,20,47,33,22,6,44,22,26,34,12,12）,

+Y=c,,,,,,,,,,,,）

>XX<-cor（cement[1:

4]）

>kappa（XX,exact=T）

[1]

>eigen（XX）

发现变量的多重共线性很强，且有

+++=0

说明X1,X2，X3，X4多重共线。

其实逐步回归可以解决多重共线的问题。

我们可以检验一下step函数去掉变量后的共线性。

step去掉了X3和X4。

我们看看去掉他们的共线性如何。

>XX<-cor（cement[1:

2]）

>kappa（XX,exact=T）

[1]

我们发现去掉X3和X4后，条件数降低好多好多。

说明step函数是合理的。

首先得把这个表格看懂。

里面的数字应该是有感染和无感染的人数。

而影响变量有三个。

我们把这些影响变量进行编码。

如下。

发生

不发生

抗生素X1

危险因子X2

有无计划X3

是否感染Y

对数据的处理，如下

频数

然后用R处理并求解模型

>hospital<（X1=rep（c（2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3）

展开阅读全文