[R] problem with ML estimation

From: souvik banerjee <bansouvik_at_gmail.com>
Date: Sun, 24 Feb 2008 16:41:55 +0530


dear list,
as a part my problem. I have to estimate some parameters using ML estimation. The form of the likelihood function is not straight forward and I had to use a for loop to define the function. I used "optim" to maximise the result but was not sure of the programme.
To validate my results, I tried to write a function to obtain the MLE of a bivariate normal in the same manner.
On generating a simulated data and running the script, I am getting results which are very inaccurate.
I am sending my programme. It would be of great help if someone can point me, where I am going wrong.

############################################################

library(MASS)
Sigma <- matrix(c(4,2.8,2.8,4),2,2)
y<-mvrnorm(n=1000, rep(0, 2), Sigma)

mlebinorm<- function(startval, y)# startval = Initial Values to be passed , y= data
{
 startval<-as.numeric(startval)
 oneside=matrix(c(0,0,0,0,-1,0,0,0,0,1),5,2)  otherside = c(-0.9999999,-0.999999)
 n<- ncol(y)
 nf<- function(x)
 {
  mu1<-x[1]
  mu2<-x[2]
  sig1<-x[3]
  sig2<-x[4]
  rho<-x[5]
  nf<-(-n/2)*(log(sig1*sig2*(1-rho^2)))-(0.5/(1-rho^2)*
(((y[1,1]-mu1)^2/sig1)-
(2*rho*(y[1,1]-mu1)*(y[1,2]-mu2))/sqrt(sig1*sig2)+
(y[1,2]-mu2)^2/sig2))

  for(i in 2:n)
  {
   nf<- nf - (0.5*(1-rho^2)*
(((y[i,1]-mu1)^2/sig1)-
(2*rho*(y[i,1]-mu1)*(y[i,2]-mu2))/sqrt(sig1*sig2)+
(y[i,2]-mu2)^2/sig2))

  }
 return(nf)

 }

 constrOptim(startval,nf,NULL,ui=t(oneside), ci=otherside,control=list(fnscale=-1))

}

mlebinorm(c(1,1,2,2,0.3),y)

############################################################

Souvik Banerjee
Lecturer
Department of statistics
Memari College
Burdwan
India

        [[alternative HTML version deleted]]



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 24 Feb 2008 - 11:16:43 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 24 Feb 2008 - 13:30:16 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive