统计建模与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] 7 10 13
> z1<-crossprod(x,y);z1
[,1]
[1,] 32
> z2<-outer(x,y);z2
[,1] [,2] [,3]
[1,] 4 5 6
[2,] 8 10 12
[3,] 12 15 18
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]
> x<-c(rep(1,5),rep(2,3),rep(3,4),rep(4,2));x
2.4
> H<-matrix(nrow=5,ncol=5)
> for (i in 1:5)
+ for(j in 1: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('要求输入一个正整数')
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-sample Kolmogorov-Smirnov test
data: x
D = 0.073, p-value = 0.6611
alternative hypothesis: two-sided
Warning message:
In ks.test(x, "pnorm", 73.668, 3.9389) :
ties should not be present for the Kolmogorov-Smirnov test
这里出现警告信息是因为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='blu e')
> 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）先求矩估计。 总体的期望为11
12
()a a a x dx a +∞
+-∞++=
+?。因此我们有
12()a E x a +=+。可解得a=(2*E （x ）-1）/(1-E(x )).因此我们用样本的均值来估计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）采用极大似然估计 首先求出极大似然函数为
L a;x = a +1 x i a =(a +1)n x i a n
i =1
n
i =1
再取对数为
ln L a;x =nln a +1 +aln( x i n
i =1
最后求导
elnL(a;x)ea =n
a +1
+ln x i n
i =1 好了下面开始用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
用极大似然估计得出λ=n/ x i n i =1.现用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)
One Sample t-test
data: x
t = -2.4534, df = 9, p-value = 0.9817
alternative hypothesis: true mean is greater than 72
95 percent confidence interval:
63.96295 Inf
sample estimates:
mean of x
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: 220 188 162 230 145 160 238 188 247 113
11: 126 245 164 231 256 183 190 158 224 175
21:
Read 20 items
> 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
这是检验两个总体是否存在差异的问题。可用符号检验和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))) (-2,0] (0,2] (2,4] (4,6] (6,8] 4 4 6 4 1 发现A中有频数小于5，故应该重新调整分组 > A<-table(cut(x,br=c(-2,2,4,8))) > A (-2,2] (2,4] (4,8] 8 6 5 然后再计算理论分布 > 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] 5 5 5 5 > 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.0 33.0 37.5 48.0 42.5 40.0 42.0 36.0 11.3 22.0 11: 36.0 27.3 14.2 32.1 52.0 38.0 17.3 20.0 21.0 46.1 21: Read 20 items > y<-scan() 1: 37.0 41.0 23.4 17.0 31.5 40.0 31.0 36.0 5.7 11.5 11: 21.0 6.1 26.5 21.3 44.5 28.0 22.6 20.0 11.0 22.3 21: Read 20 items > 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) 结果显示这两种疗法没什么区别