[求助]SPLUS中极大似然估计问题！

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我在用SPLUS做极大似然估计的过程中，已经得到似然函数L（p,q,r,s,t),下一步是分别对p,q,r,s,t求偏导数取零来联立一个方程组，从而解得5个参数得估计值，我不知道该如何在软件中用S语言实现这一步骤，急！希望高手指点！不胜感激！！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：极大似然估计 极大似然 PLUS 似然估计 Plu splus 极大似然估计

急啊～～

那位大哥帮帮忙！

S-plus没用过，但是和R都是用S语言，所以R的我贴出来，看看对大家有没有帮助...


1. 数据与模型

我们要使用的数据来自于“MASS”包中的geyser数据。先把数据调出来，看看它长什么样子。

> geyser
    waiting  duration
1        80 4.0166667
2        71 2.1500000
3        57 4.0000000
4        80 4.0000000
5        75 4.0000000
......
该数据采集自美国黄石公园内的一个名叫Old Faithful 的喷泉。“waiting”就是喷泉两次喷发的间隔时间，“duration”当然就是指每次喷发的持续时间。在这里，我们只用到“waiting”数据，为了简单一点，可以使用attach()函数。

> attach(geyser)

2. 模型

绘制出数据的频率分布直方图：

> hist(waiting)
ml_hist


从图中可以看出，其分布是两个正态分布的混合。可以用如下的分布函数来描述该数据

f(x)=pN(xi;μ1,σ1)+(1−p)N(xi;μ2,σ2)
该函数中有5个参数[latex p、\mu_1、\sigma_1、\mu_2、\sigma_2$需要确定。上述分布函数的对数极大似然函数为：

l=∑ni=1log{pN(xi;μ1,σ1)+(1−p)N(xi;μ2,σ2)}
   
3. 估计

3.1 在R中定义对数似然函数：

> #定义log-likelihood函数
> LL<-function(params,data)
+ {#参数"params"是一个向量，依次包含了五个参数：p,mu1,sigma1,
+ #mu2,sigma2.
+ #参数"data"，是观测数据。
+ t1<-dnorm(data,params[2],params[3])
+ t2<-dnorm(data,params[4],params[5])
+ #这里的dnorm()函数是用来生成正态密度函数的。
+ f<-params[1]*t1+(1-params[1])*t2
+ #混合密度函数
+ ll<-sum(log(f))
+ #log-likelihood函数
+ return(-ll)
+ #nlminb()函数是最小化一个函数的值，但我们是要最大化log-
+ #likeilhood函数，所以需要在“ll”前加个“-”号。
+ }
  
3.2 参数估计

> #用hist函数找出初始值
> hist(waiting,freq=F)
> lines(density(waiting))
> #拟合函数####optim####
> geyser.res<-nlminb(c(0.5,50,10,80,10),LL,data=waiting,
+ lower=c(0.0001,-Inf,0.0001,-Inf,-Inf,0.0001),
+ upper=c(0.9999,Inf,Inf,Inf,Inf))
> #初始值为p=0.5,mu1=50,sigma1=10,mu2=80,sigma2=10
> #LL是被最小化的函数。
> #data是拟合用的数据
> #lower和upper分别指定参数的上界和下界。
  
3.3 估计结果

> #查看拟合的参数
> geyser.res$par
[1] 0.3075937 54.2026518 4.9520026 80.3603085 7.5076330
> #拟合的效果
> X<-seq(40,120,length=100)
> #读出估计的参数
> p<-geyser.res$par[1]
> mu1<-geyser.res$par[2]
> sig1<-geyser.res$par[3]
> mu2<-geyser.res$par[4]
> sig2<-geyser.res$par[5]
> #将估计的参数函数代入原密度函数。
> f<-p*dnorm(X,mu1,sig1)+(1-p)*dnorm(X,mu2,sig2)
> #作出数据的直方图
> hist(waiting,probability=T,col=0,ylab="Density",
+ ylim=c(0,0.04),xlab="Eruption waiting times")
> #画出拟合的曲线
> lines(X,f)
clip_image004.jpg

> detach()
  
小结

在R中作极大似然估计，主要就是定义似然后函数，然后再用nlminb函数对参数进行估计。

参考文献：

l Brian S. Everitt(2002). A Handbook of Statistical Analyses Using S-Plus(Second Edition). CRC Press LLC