R语言实验6.docx
《R语言实验6.docx》由会员分享，可在线阅读，更多相关《R语言实验6.docx（15页珍藏版）》请在冰点文库上搜索。
R语言实验6
R语言实验6
R语言实验6
令
ln
=0
解此似然方程得到
，或写为
。
容易验证
，从而α使得L达到极大，即参数α的极大似然估计量un
。
以下请根据上式完成R程序，计算出参数α的极大似然估计量
的值。
源代码及运行结果：
（复制到此处，不需要截图）
>f<-function（a）6/（a+1）+sum（log（x））
>uniroot（f,c（0,1））
$root
[1]0.211182
$f.root
[1]-3.844668e-05
$iter
[1]5
$init.it
[1]NA
$estim.prec
[1]6.103516e-05
1.（习题4.2）设元件无故障工作时间X具有指数分布，取1000个元件工作时间的记录数据，经分组后得到它的频数分布为
组中值xi
5
15
25
35
45
55
65
频数vi
365
245
150
100
70
45
25
如果各组中数据都取为组中值，试用极大似然函数估计求λ的点估计。
提示：
①根据教材P168例4.7知，指数分布中参数λ的极大似然估计是n/
。
②利用rep（）函数。
解：
源代码及运行结果：
（复制到此处，不需要截图）
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）
[1]0.05
2.（习题4.3）为检验某自来水消毒设备的效果，现从消毒后的水中随机抽取50升，化验每升水中大肠杆菌的个数（假设一升水中大肠杆菌个数服从Poisson分布），其化验结果如下：
大肠杆菌数/升
0
1
2
3
4
5
6
水的升数
17
20
10
2
1
0
0
试问平均每升水中大肠杆菌个数为多少时，才能使上述情况的概率为最大？
解：
此题实际就是求泊松分布中参数λ的极大似然估计。
泊松分布的分布律为P{X=k}=
k=0,1,2,…,λ>0
设x1，x2，…，xn为其样本X1，X2，…，Xn的一组观测值。
于是此分布的似然函数为L（λ;x）=L（λ;x1，…，xn）=
相应的对数似然函数为lnL（λ;x）=-nλ+
lnλ-
令
=0
解此似然方程得到
容易验证
，从而λ使得L达到极大，即参数λ的极大似然估计量
。
以下请据此完成R程序，计算出参数λ的极大似然估计量
的值。
同上题，也需要利用rep（）函数。
源代码及运行结果：
（复制到此处，不需要截图）
x<-c（rep（0,17）,rep（1,20）,rep（2,10）,rep（3,2）,rep（4,1）,rep（5,0）,rep（6,0））
>mean（x）
[1]1
3.（习题4.4）利用R软件中的nlm（）函数求解无约束优化问题
minf（x）=（-13+x1+（（5-x2）x2-2）x2）2+（-29+x1+（（x2+1）x2-14）x2）2,
取初始点x（0）=（0.5,-2）T。
提示：
参考教材P173公式（4.13）对应的例题。
解：
源代码及运行结果：
（复制到此处，不需要截图）
obj<-function（x）{
f<-c（-13+x[1]+（（5-x[2]）*x[2]-2）*x[2],-29+x[1]+（（x[2]+1）*x[2]-14）*x[2]）
sum（f^2）
}
>source（"Rosenbrouk.R"）
>x0<-c（0.5,-2）
>nlm（obj,x0）
$minimum
[1]48.98425
$estimate
[1]11.4127791-0.8968052
$gradient
[1]1.413268e-08-1.462297e-07
$code
[1]1
$iterations
[1]16
结论：
$minimum是函数的最目标值，即f*=48.98425，$estimate是最优点的估计值，即x*=（11.4127791,-0.8968052）T;$gradient是在最优点处（估计值）目标函数梯度值，即f*（1.413268e-08,-1.462297e-07）T;$code是指标，这里是1，表示迭代成功；$iterations
是迭代次数，这里是16，表示迭代了16次迭代。
4.（习题4.5）正常人的脉搏平均每分钟72次，某医生测得10例四乙基铅中毒患者的脉搏数（次/min）如下：
54676878706667706569
已知人的脉搏次数服从正态分布，试计算这10名患者平均脉搏次数的点估计和95%的区间估计。
并作单侧区间估计，试分析这10患者的平均脉搏次数是否低于正常人的平均脉搏次数。
提示：
此题是一个正态总体的估计问题，且由于总体方差未知，因此可以直接使用R语言中t.test（）函数进行分析。
解：
源代码及运行结果：
（复制到此处，不需要截图）
x<-c（54,67,68,78,70,66,67,70,65,69）
>t.test（x）#t.rest（）做单样本正态分布区间估计
OneSamplet-test
data:
x
t=35.947,df=9,p-value=4.938e-11
alternativehypothesis:
truemeanisnotequalto0
95percentconfidenceinterval:
63.158571.6415
sampleestimates:
meanofx
67.4
>t.test（x,alternative="less"，mu=72）#t.rest（）做单样本正态分布单侧间估计
OneSamplet-test
data:
x
t=-2.4534,df=9,p-value=0.01828
alternativehypothesis:
truemeanislessthan72
95percentconfidenceinterval:
-Inf70.83705
sampleestimates:
meanofx
67.4
结论：
双侧区间估计：
平均脉搏点估计为67.4，95%区间估计为63.158571.6415；
单侧区间估计：
p=0.01828<0.05,拒绝原假设，平均脉搏低于正常人。
5.（习题4.6）甲、乙两种稻种分别播种在10块试验田中，每块试验田甲、乙稻种各种一半。
假设两稻种产量X,Y均服从正态分布，且方差相等。
收获后10块试验田的产量如下所示（单位：
千克）。
甲种
140
137
136
140
145
148
140
135
144
141
乙种
135
118
115
140
128
131
130
115
131
125
求出两稻种产量的期望差μ1-μ2的置信区间（α=0.05）。
提示：
此题是两个正态总体的区间估计问题，且由于两总体方差未知，因此可以直接使用R语言中t.test（）函数进行分析。
t.test（）可做两正态样本均值差的估计。
注意此例中两样本方差相等。
参考教材P185倒数第四行开始至P186。
解：
源代码及运行结果：
（复制到此处，不需要截图）
>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=TRUE）
TwoSamplet-test
data:
xandy
t=4.6287,df=18,p-value=0.0002087
alternativehypothesis:
truedifferenceinmeansisnotequalto0
95percentconfidenceinterval:
7.53626120.063739
sampleestimates:
meanofxmeanofy
140.6126.8
结论：
期望差的95%置信区间为7.53626120.063739
6.（习题4.7）甲、乙两组生产同种导线，现从甲组生产的导线中随机抽取4根，从乙组生产的导线中随机抽取5根，它们的电阻值（单位：
Ω）分别为
甲组
0.143
0.142
0.143
0.137
已组
0.140
0.142
0.136
0.138
0.140
假设两组电阻值分别服从正态分布N（μ1,σ2）和N（μ1,σ2），σ2未知。
试求μ1-μ2的置信区间系数为0.95的区间估计。
提示：
此题是两个正态总体的估计问题，且由于两总体方差未知，因此可以直接使用R语言中t.test（）函数进行分析。
t.test（）可做两正态样本均值差的估计。
注意此例中两样本方差相等。
参考教材P185倒数第四行开始至P186。
解：
源代码及运行结果：
（复制到此处，不需要截图）
>x<-c（0.143,0.142,0.143,0.137）
>y<-c（0.140,0.142,0.136,0.138,0.140）
>t.test（x,y,var.equal=TRUE）
TwoSamplet-test
data:
xandy
t=1.198,df=7,p-value=0.2699
alternativehypothesis:
truedifferenceinmeansisnotequalto0
95percentconfidenceinterval:
-0.0019963510.006096351
sampleestimates:
meanofxmeanofy
0.141250.13920
结论：
期望差的95%区间估计为-0.0019963510.006096351
7.（习题4.8）对习题4.6中甲乙两种稻种的数据作方差比的区间估计，并用其估计值来判定两数据是否等方差。
若两数据方差不相等，试重新计算两稻种产量的期望差μ1-μ2的置信区间（α=0.05）。
提示：
在R软件中，var.test（）函数能够提供两个样本方差比的区间估计。
参考教材P189倒数第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）
>var.test（x,y）
Ftesttocomparetwovariances
data:
xandy
F=0.23533,numdf=9,denomdf=9,p-value=0.04229
alternativehypothesis:
trueratioofvariancesisnotequalto1
95percentconfidenceinterval:
0.058452760.94743902
sampleestimates:
ratioofvariances
0.2353305
>t.test（x,y,）
WelchTwoSamplet-test
data:
xandy
t=4.6287,df=13.014,p-value=0.0004712
alternativehypothesis:
truedifferenceinmeansisnotequalto0
95percentconfidenceinterval:
7.35971320.240287
sampleestimates:
meanofxmeanofy
140.6126.8
结论：
期望差的95%置信区间为7.35971320.240287
8.（习题4.9）设电话总机在某段时间内接到的呼唤的次数服从参数未知的Poisson分布P（λ），现收集了42个数据
接到呼唤次数
0
1
2
3
4
5
6
出现的频数
7
10
12
8
3
2
0
试求出平均呼唤次数λ的估计值和它的置信系数为0.95的置信区间。
此题给出完整答案，供参考。
解：
从习题4.3已经知道，平均呼唤次数λ的估计值就是⎺X。
此题是非正态总体的区间估计问题。
但由中心极限定理可知，此类问题其置信区间与方差σ2未知时一个正态总体的均值的置信区间完全一样，即为
。
参考教材P190公式（4.46）及其后的程序。
源代码及运行结果：
>x<-c（rep（0,7）,rep（1,10）,rep（2,12）,rep（3,8）,rep（4,3）,rep（5,2））
>n<-length（x）
>tmp<-sd（x）/sqrt（n）*qnorm（1-0.05/2）
>mean（x）
[1]1.904762
>mean（x）-tmp;mean（x）+tmp#置信区间的下限和上限
[1]1.494041
[1]2.315483
结论：
平均呼唤次数为1.9
0.95的置信区间为[1.49,2.32]。
9.（习题4.10）已知某种灯泡寿命服从正态分布，在某星期所生产的该灯泡中随机抽取10只，测得其寿命（单位：
小时）为
1067919119678511269369181156920948
求灯泡寿命平均值的置信度为0.95的单侧置信下限。
提示：
此题是一个正态总体的区间估计问题，且由于总体方差未知，因此可以直接使用R语言中t.test（）函数进行分析。
根据教材P191定义4.7知，单侧置信下限，t.test（）函数中的参数alternative="greater"。
参考教材P193-194的例4.22及其后面一段文字。
解：
源代码及运行结果：
（复制到此处，不需要截图）
x<-c（1067,919,1196,785,1126,936,918,1156,920,948）
>t.test（x,alternative="greater"）
OneSamplet-test
data:
x
t=23.969,df=9,p-value=9.148e-10
alternativehypothesis:
truemeanisgreaterthan0
95percentconfidenceinterval:
920.8443Inf
sampleestimates:
meanofx
997.1
结论：
灯泡平均寿命置信度95%的单侧置信度下限为920.8443
思考：
1.常用的点估计的方法有哪些？
距法，极大似然法，最小二乘法等。
2.估计量的优良性准则有哪些？
无偏性，有效性和相合性（一致性）
3.在R语言中，求矩法估计时，如果令总体的k阶原点矩等于它样本的k阶原点矩得到的k个方程无法求出解析解，此时需要调用教材P103的2.9.3节作者编写的Newton法程序求此方程组的解。
在求极大似然估计时需要求解（对数）似然函数的极大值，如果无法求出其相应似然方程的解析解，若此（对数）似然函数是一维变量，可以用R中的_optimeize（）（或optimise（））_____函数变相求解此（对数）似然函数的相反数的极小值？
若此（对数）似然函数是多维变量，可以用R中的__nlm（）__函数变相求解
–L（θ;x）或
–lnL（θ;x）？
一、实验小结（必写，但字数不限）
这次实验主要学会掌握矩法估计与极大似然估计的求法，会利用R软件完成一个和两个正态总体的区间估计，非正态总体的区间估计和单侧置信区间估计，还要会分析结果，得出结论。
懂得什么样的题，用什么样的方法和函数。
通过这次实验，已经掌握了一些基本参数估计知识，要多看书本，多加练习。
升级会员 