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]]

