最新范文R语言 实验6 参数估计.docx
《最新范文R语言 实验6 参数估计.docx》由会员分享,可在线阅读,更多相关《最新范文R语言 实验6 参数估计.docx(14页珍藏版)》请在冰豆网上搜索。
最新范文R语言实验6参数估计
R语言实验6参数估计
一、实验目的:
1.掌握矩法估计与极大似然估计的求法;
2.学会利用R软件完成一个和两个正态总体的区间估计;
3.学会利用R软件完成非正态总体的区间估计;
4.学会利用R软件进行单侧置信区间估计。
二、实验内容:
练习:
要求:
①完成练习并粘贴运行截图到文档相应位置(截图方法见下),并将所有自己输入文字的字体颜色设为红色(包括后面的思考及小结),②回答思考题,③简要书写实验小结。
④修改本文档名为“本人完整学号姓名1”,其中1表示第1次实验,以后更改为2,3,...。
如文件名为“1305543109张立1”,表示学号为1305543109的张立同学的第1次实验,注意文件名中没有空格及任何其它字符。
最后连同数据文件、源程序文件等(如果有的话,本次实验没有),一起压缩打包发给课代表,压缩包的文件名同上。
截图方法:
法1:
调整需要截图的窗口至合适的大小,并使该窗口为当前激活窗口(即该窗口在屏幕最前方),按住键盘Alt键(空格键两侧各有一个)不放,再按键盘右上角的截图键(通常印有“印屏幕”或“PrScrn”等字符),即完成截图。
再粘贴到word文档的相应位置即可。
法2:
利用QQ输入法的截屏工具。
点击QQ输入法工具条最右边的“扳手”图标
选择
其中的“截屏”工具。
)
1.自行完成教材P163页开始的4.1.3-4.3节中的例题。
2.(习题4.1)设总体的分布密度函数为⎩⎨⎧<<+=,
0,10)1();(其他xxxfα
ααX1,X2,…,Xn为其样本,求参数α的矩估计量1ˆα
和极大似然估计量2ˆα。
现测得样本观测值为
0.1,0.2,0.9,0.8,0.7,0.7
求参数α的估计值。
解:
先求参数α的矩估计量1ˆα
。
由于只有一个参数,因此只需要考虑E(X)=X。
而由E(X)的定义有:
E(X)=
2
1|21)1()(1021
0++=++=+⋅=⋅++∞∞-⎰⎰αααααααxdxxxdxxfx因此X=++21αα,解得211ˆ1--=X
α。
以下请根据上式完成R程序,计算出参数α的矩估计量1ˆα的值。
源代码及运行结果:
(复制到此处,不需要截图)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ˆα
。
只需要考虑x∈(0,1)部分。
依题意,此分布的似然函数为L(α;x)=
∏∏==+=n
ii
n
nii
xxf1
1
)()1();(ααα
相应的对数似然函数为lnL(α;x)=nln(α+1)+αln
∏=n
ii
x
1
令
++=∂∂1);(lnαααn
xLln∏=n
iix1
=0解此似然方程得到1ln1--
=∏=n
ii
xnα,或写为1ln1
--
=∑=n
ii
x
n
α。
容易验证0ln2
2<∂∂αL
从而α使得L达到极大,即参数α的极大似然估计量un1lnˆ1
2--=∑=n
ii
X
n
α
。
以下请根据上式完成R程序,计算出参数α的极大似然估计量2ˆα
的值。
源代码及运行结果:
(复制到此处,不需要截图)>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
3.(习题
4.2)设元件无故障工作时间X具有指数分布,取1000个元件工作时间的
提示:
①根据教材P168例4.7知,指数分布中参数λ的极大似然估计是n/
∑=n
ii
X
1
。
②利用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
4.(习题4.3)为检验某自来水消毒设备的效果,现从消毒后的水中随机抽取50升,
化验每升水中大肠杆菌的个数(假设一升水中大肠杆菌个数服从
Poisson分布),其化验结果如下:
试问平均每升水中大肠杆菌个数为多少时,才能使上述情况的概率为最大?
解:
此题实际就是求泊松分布中参数λ的极大似然估计。
泊松分布的分布律为P{X=k}=
!
kekλ
λ-,k=0,1,2,…,λ>0
设x1,x2,…,xn为其样本X1,X2,…,Xn的一组观测值。
于是此分布的似然函数为L(λ;x)=L(λ;x1,…,xn)=
λλ
λ
λnnxn
iixexxxen
ii
i
-=-∑=
=∏
!
!
!
11
1
相应的对数似然函数为lnL(λ;x)=-nλ+
∑=n
ii
x
1
lnλ-
∑=n
ii
x1
)!
ln(
令∑=++-=∂∂n
iixnxL1
1);(lnλαλ=0
解此似然方程得到x=λ
容易验证0ln2
2<∂∂λ
L,从而λ使得L达到极大,即参数λ的极大似然估计量X=λˆ。
以下请据此完成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
5.(习题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次迭代。
6.(习题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,拒绝原假设,平均脉搏低于正常人。
7.(习题4.6)甲、乙两种稻种分别播种在10块试验田中,每块试验田甲、乙稻种各
种一半。
假设两稻种产量X,Y均服从正态分布,且方差相等。
收获后10块试验田
12
提示:
此题是两个正态总体的区间估计问题,且由于两总体方差未知,因此可以直接使用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
8.(习题4.7)甲、乙两组生产同种导线,现从甲组生产的导线中随机抽取4根,从
)分别为
假设两组电阻值分别服从正态分布N(1,)和N(1,),未知。
试求μ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
9.(习题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
10.(习题4.9)设电话总机在某段时间内接到的呼唤的次数服从参数未知的Poisson
此题给出完整答案,供参考。
解:
从习题4.3已经知道,平均呼唤次数λ的估计值就是⎺X。
此题是非正态总体的区间估计问题。
但由中心极限定理可知,此类问题其置信区间与方差σ2未知时一个正态总体的均值的置信区间完全一样,即为]+,[2/2/ααZnS
XZnS
X-。
参考教材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]。
11.(习题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()__函数变相求解θmin–L(θ;x)或θ
min–lnL(θ;x)?
三、实验小结(必写,但字数不限)
这次实验主要学会掌握矩法估计与极大似然估计的求法,会利用R软件完成一个和两个正态总体的区间估计,非正态总体的区间估计和单侧置信区间估计,还要会分析结果,得出结论。
懂得什么样的题,用什么样的方法和函数。
通过这次实验,已经掌握了一些基本参数估计知识,要多看书本,多加练习。