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