当前位置:文档之家› R语言实验6

R语言实验6

实验6 参数估计一、实验目的:1. 掌握矩法估计与极大似然估计的求法;2. 学会利用R 软件完成一个和两个正态总体的区间估计;3. 学会利用R 软件完成非正态总体的区间估计;4. 学会利用R 软件进行单侧置信区间估计。

二、实验内容:练习: 要求:①完成练习并粘贴运行截图到文档相应位置(截图方法见下),并将所有自己输入文字的字体颜色设为红色(包括后面的思考及小结),②回答思考题,③简要书写实验小结。

④修改本文档名为“本人完整学号姓名1”,其中1表示第1次实验,以后更改为2,3,...。

如文件名为“1305543109张立1”,表示学号为1305543109的张立同学的第1次实验,注意文件名中没有空格及任何其它字符。

最后连同数据文件、源程序文件等(如果有的话,本次实验没有),一起压缩打包发给课代表,压缩包的文件名同上。

截图方法:法1:调整需要截图的窗口至合适的大小,并使该窗口为当前激活窗口(即该窗口在屏幕最前方),按住键盘Alt 键(空格键两侧各有一个)不放,再按键盘右上角的截图键(通常印有“印屏幕”或“Pr Scrn ”等字符),即完成截图。

再粘贴到word 文档的相应位置即可。

法2:利用QQ 输入法的截屏工具。

点击QQ 输入法工具条最右边的“扳手”图标,选择其中的“截屏”工具。

)1. 自行完成教材P163页开始的4.1.3-4.3节中的例题。

2. (习题4.1)设总体的分布密度函数为⎩⎨⎧<<+=,,10)1();(其他x x x f αααX 1,X 2,…,X n 为其样本,求参数α 的矩估计量1ˆα和极大似然估计量2ˆα。

现测得样本观测值为0.1, 0.2, 0.9, 0.8, 0.7, 0.7求参数 α 的估计值。

解:先求参数α 的矩估计量1ˆα。

由于只有一个参数,因此只需要考虑E(X )=X 。

而由E(X )的定义有:E(X )=21|21)1()(1021++=++=+⋅=⋅++∞∞-⎰⎰αααααααx dx x x dx x f x 因此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 )=∏∏==+=ni inn i ix x f 11)()1();(ααα相应的对数似然函数为 ln L (α; x )=n ln(α +1)+ α ln∏=ni ix1令 ++=∂∂1);(ln αααnx L ln ∏=ni i x 1=0解此似然方程得到1ln 1--=∏=ni ix n α,或写为1ln 1--=∑=ni ixnα。

容易验证0ln 22<∂∂αL,从而α 使得L 达到极大,即参数α 的极大似然估计量un 1ln ˆ12--=∑=ni iXnα。

以下请根据上式完成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-053. (习题4.2)设元件无故障工作时间X 具有指数分布,取1000个元件工作时间的提示:①根据教材P168例4.7知,指数分布中参数 λ 的极大似然估计是n /∑=ni iX1。

②利用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.054. (习题4.3)为检验某自来水消毒设备的效果,现从消毒后的水中随机抽取50升,化验 每升水中大肠杆菌的个数(假设一升水中大肠杆菌个数服从Poisson 分布),其化验结果如下:试问平均每升水中大肠杆菌个数为多少时,才能使上述情况的概率为最大? 解:此题实际就是求泊松分布中参数 λ 的极大似然估计。

泊松分布的分布律为 P {X =k }=!k e k λλ-, k =0,1,2,…, λ > 0设x 1,x 2,…,x n 为其样本X 1,X 2,…,X n 的一组观测值。

于是此分布的似然函数为 L (λ; x )= L (λ; x 1,…,x n )=λλλλn n x ni i x e x x x e ni ii-=-∑==∏!!!111相应的对数似然函数为 ln L (λ; x )= -n λ+∑=ni ix1ln λ-∑=ni ix 1)!ln(令∑=++-=∂∂ni i x n x L 11);(ln λαλ=0 解此似然方程得到x =λ容易验证0ln 22<∂∂λ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] 15.(习题4.4)利用R软件中的nlm()函数求解无约束优化问题min f(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)如下:54 67 68 78 70 66 67 70 65 69已知人的脉搏次数服从正态分布,试计算这10名患者平均脉搏次数的点估计和95%的区间估计。

并作单侧区间估计,试分析这10患者的平均脉搏次数是否低于正常人的平均脉搏次数。

提示:此题是一个正态总体的估计问题,且由于总体方差未知,因此可以直接使用R 语言中t.test()函数进行分析。

解:源代码及运行结果:(复制到此处,不需要截图)x<-c(54 , 67 , 68 , 78 ,70 , 66 , 67 , 70 , 65, 69)> t.test(x) #t.rest()做单样本正态分布区间估计One Sample t-testdata: xt = 35.947, df = 9, p-value = 4.938e-11alternative hypothesis: true mean is not equal to 095 percent confidence interval:63.1585 71.6415sample estimates:mean of x67.4> t.test(x,alternative="less",mu=72) #t.rest()做单样本正态分布单侧间估计One Sample t-testdata: xt = -2.4534, df = 9, p-value = 0.01828alternative hypothesis: true mean is less than 7295 percent confidence interval:-Inf 70.83705sample estimates:mean of x67.4结论:双侧区间估计:平均脉搏点估计为67.4,95%区间估计为63.1585 71.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)Two Sample t-testdata: x and yt = 4.6287, df = 18, p-value = 0.0002087alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:7.536261 20.063739sample estimates:mean of x mean of y140.6 126.8结论:期望差的95%置信区间为7.536261 20.0637398.(习题4.7)甲、乙两组生产同种导线,现从甲组生产的导线中随机抽取4根,从)分别为假设两组电阻值分别服从正态分布N(1, )和N(1, ),未知。

相关主题