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@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.