2
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ Please attach the profile of the likelihood from 1 to the unconstrained optimum. Does "sensible" mean the second parameter is between 0 and 1? Then the algorithm just needs to explore outside the interval for a few steps but the restriction is unnecessary. $\endgroup$
    – Alex
    Jul 25, 2011 at 15:05
  • 1
    $\begingroup$ your code prints out a bug 'loglik.negbin.perageclass(starting.values)' yields 'Error: object 'C_a' not found'. Correct please. $\endgroup$
    – user603
    Jul 25, 2011 at 16:30

0

Browse other questions tagged or ask your own question.