calculate_mle<-function(x,n0,eta,v,sigma2){
log_likelihoodFun<-function(theta){
eta<-theta[1]
v<-theta[2]
sigma2<-theta[3]
delta<-1-exp(v)
if(sigma2 < 0 | is.nan(sigma2) | delta< 0 | delta >1 ) return(NA)
logL<-n0*log(1-exp(v))+length(x)*v-(length(x)/2)*log(2*pi)-(length(x)/2)*log(sigma2)-((1/(2*sigma2))*sum((x-eta+v+0.5*sigma2)*(x-eta+v+0.5*sigma2)))
return(-logL)
}
fit<-optim(par = c(eta,v,sigma2),fn = log_likelihoodFun)
return(list(etahat=fit$par[1],vhat=fit$par[2],sigma2hat=fit$par[3]))
}
delta<-0.1
sigma2<-0.1
mu<-0
n<-30
n0<-0
n1<-n-n0
x<-rnorm(n1,mu,sqrt(sigma2))
x_bar<-mean(x)
ss<-var(x)
delta_mle<-n0/n
sigma2_mle<-((n1-1)/n1)*ss
mle<-calculate_mle(x,n0,log(1-delta_mle)+x_bar+0.5*sigma2_mle,log(1-delta_mle),sigma2_mle)
相似问题