For my research I am currently working in R and I have created a function which gives me the loglikelihood. Now, I would like to optimize it for multiple parameters. So I decided to use the optim() function in R.
Since I have 1 vector parameter and another parameter which are constrained to have values between 0 and 1, since they are proportions, I used the option method="L-BFGS-B".
Applying this I get results but one of my restricted parameters is estimated as exactly 1.0000. If I do not constrain it, it gives a sensible answer. But, the standard error is large (still smaller than under the constrained). Is there a way to remain with the restriction but that I can solve the fact that the program ends with 1.000?
I give some piece of my code.
Ps<-1-c(0.02083333, 0.03067485, 0.11290323, 0.13684211, 0.23629490)
loglik.negbin.perageclass<-function(param.qphi){
#q.param<-rep(param.qphi[1],5)
#phi.param<-rep(param.qphi[2],5)
#size.param<-param.qphi[3]
q.param<-param.qphi[1:5]
phi.param<-rep(param.qphi[6],5)
size.param<-param.qphi[7]
shape.param<-param.qphi[8]
scale.param<-param.qphi[9]
#I_days[1,]<-param.qphi[8:12]
delta.touse<-dgamma(0:4,shape.param,scale.param)/sum(dgamma(0:4,shape.param,scale.param))
delta<-function(t_current){
diff<-rep(0,t_current)
res<-matrix(0,nrow=t_current,ncol=5)
for(k in 1:t_current){
diff[k]=t_current-k
if(diff[k]==0){
res[k,]<-delta.touse[1]
}
else if(diff[k]==1){
res[k,]<-delta.touse[2]
}
else if(diff[k]==2){
res[k,]<-delta.touse[3]
}
else if(diff[k]==3){
res[k,]<-delta.touse[4]
}
else if(diff[k]==4){
res[k,]<-delta.touse[5]
}
}
return(res)
}
week<-1
for(k in 2:35){
I_a_days[k,]<-phi.param*(q.param*C_a[,])%*%(apply(delta(k-1)*Ps*N*I_days[1:k-1,],2,sum))
I_s_days[k,]<-(1-phi.param)*(q.param*C_s[,])%*%(apply(delta(k-1)*Ps*N*I_days[1:k-1,],2,sum))
I_days[k,]<-I_a_days[k,]+I_s_days[k,]
if(k%in%newweek){
beginningweek<-k-6
for(l in beginningweek:k){
I_week[week,]<-I_week[week,]+I_days[l,]
I_s_week[week,]<-I_s_week[week,]+I_s_days[l,]
I_a_week[week,]<-I_a_week[week,]+I_a_days[l,]
}
week<-week+1
}
}
print(I_week)
# Plotting
plot(GP.data$Weekno[24:29],log(GP.data$X0.4[24:29]),pch=19,lwd=2,main="GP data",ylim=range(0,log(700)),col=1,xlab="week",ylab="number of cases")
points(GP.data$Weekno[24:29],log(GP.data$X5.14[24:29]),pch=19,col=2)
points(GP.data$Weekno[24:29],log(GP.data$X15.44[24:29]),pch=19,col=3)
points(GP.data$Weekno[24:29],log(GP.data$X45.64[24:29]),pch=19,col=4)
points(GP.data$Weekno[24:29],log(GP.data$X65.[24:29]),pch=19,col=5)
for (i in 1:5){lines(c(200925:200929),log(I_week[,i]),col=i,lwd=2)}
for (i in 1:5){lines(c(200925:200929),log(I_week[,i]),col=i,lwd=2)}
# End plotting
we<-5/t(t(GP.matr[25:29,13:17])/N)/apply(1/t(t(GP.matr[25:29,13:17])/N),1,sum)
ll<-rep(0,5)
mu.param<-rep(0,5)
for(a in 1:5){
gp.a<-a+4
mu.param<-I_week[,a]
#print(c("debug",N,N[a]*GP.matr[25:29,gp.a]))
ll[a]<-sum(we[,a],log(1E-7+dnbinom(GP.matr[25:29,gp.a],mu=mu.param,size=size.param)))
}
ll.report<-sum(ll)
return(-ll.report)
}
q_start<-rep(0.05,5)
phi_start<-0.10
size_start<-1.2
shape_start<-2.9
scale_start<-0.5
starting.values<-c(q_start,phi_start,size_start,shape_start,scale_start)
loglik.negbin.perageclass(starting.values)
res<-optim(starting.values,loglik.negbin.perageclass,method="L-BFGS- B",lower=c(0,0,0,0,0),upper=c(1,1,Inf,Inf,Inf),
control=list(trace=TRUE,maxit=500),hessian=T)
Does anyone has any idea if there is a solution for this?
Thank you very much! Kim