How to add constraint into loglikelihood function?

Question

I have a time series model (INGARCH):

lambda_t = alpha0 + alpha1*(x_(t-1)) + beta1*(lambda_(t-1))

X_t ~ poisson (lambda_t)

where t is the length of observation or data, alpha0, alpha1 and beta1 are the parameters.

X_t is the series of data, lambda_t is the series of mean.

This model has the condition of alpha1 + beta1 < 1.

In my estimation, I would like to add in the condition of alpha1 + beta1 <1 in my code, I add a while loop in the log-likelihood function, but the loop cannot stop.

What could I do to solve this problem? Is there any other way to add a constraint of alpha1 + beta1 < 1 without using while loop?

Below are my code:

ll <- function(par) {
  h.new  = rep(0,n)
  #par[1] is alpha0 
  #par[2] is alpha1
  #par[3] is beta1
  while(par[2] + par[3] < 1){
  for (i in 2:n) {
    h.new[i] <- par[1] + par[2] * dat[i-1] + par[3] * h.new[i-1]

  }
  -sum(dpois(dat, h.new, log=TRUE))
  }
}

#simply generate a dataset as I have not found a suitable real dataset to fit in
set.seed(77)
n=400
dat <- rpois(n,36)

nlminb(start = c(0.1,0.1,0.1), lower = 1e-6, ll)

You are welcome, but no need to post a thanks in the answer. Indicate helpfulness by up-voting and accepting the answer so other users see it is useful. Good luck with the rest of your questions. — kabanus, Nov 27 '18 at 09:19

kabanus · Accepted Answer · 2018-11-27T09:20:55.057

1

You do not change par at all inside the while. In particular, if you would have printed par[1] and par[2] in the while you would see you are endlessly printing the original values, 0.1 - hence you are stuck in the while for ever.

par is a single, unchanging object in each call from nlminb. You just have to make sure if par is bad, you return something not minimal, so nlminb does not keep searching in that direction:

ll <- function(par) {       
    #If alpha + beta > 1, this is terrible and return an infinite score
    #It may be better to throw an error if you get NaN values! The if will
    #fail anyway, but if you want to power through add checks:
    if( is.nan(par[2]) || is.nan(par[3]) || par[2]+par[3]>1) return(Inf)
    h.new  = rep(0,n)
    #remove while
    for (i in 2:n) {
        h.new[i] <- par[1] + par[2] * dat[i-1] + par[3] * h.new[i-1]
    }
    -sum(dpois(dat, h.new, log=TRUE))
}

The algorithm nlminb (or any minimization function) very roughly goes:

Set parameters to initial guess
Send parameters to the objective functions
Guess new parameters:

a. if the score did not improve much, return minimized guess

b. if the score is good, keep searching in this direction

c. else, search in some other direction
Go back to (2) with new parameters

Note you have to return a score for each set of parameters, you do not iterate them in the objective function.

edited Nov 27 '18 at 09:20

answered Nov 27 '18 at 08:43

kabanus

24,623
6
41
74

Thank you for your comment, @kabanus. is there any way to restrict par[2] + par[3] < 1 ? – Miyazaki Nov 27 '18 at 08:50
@Miyazaki The answer shows exactly how to do this. Please read carefully everything. – kabanus Nov 27 '18 at 08:51
In particular, the first line I added your function does this. – kabanus Nov 27 '18 at 08:51
thank you for your clarification, I have run the code again, but it returns Error in if (par[2] + par[3] > 1) return(inf) : missing value where TRUE/FALSE needed how to avoid the error? – Miyazaki Nov 27 '18 at 08:56
@Miyazaki you are getting `NaN` values likely, added a check. – kabanus Nov 27 '18 at 09:02
I have run the code you modified, and found that the optimization stop at 1 iteration and return $`par` [1] 0.1 0.1 0.1 $objective [1] Inf – Miyazaki Nov 27 '18 at 09:14
@Miyazaki that is because your score function does not work properly (that is why you got NaN values in the first place). Fixing that is out of scope for this question. I answered your question - how to add a constraint. If you want to ask a new question on why does your ll function not provide a good result then go ahead, but it is not polite to change the original question. If this answers that - accept and post any other new question you may have. This is the way SO works and will provide you the best results. – kabanus Nov 27 '18 at 09:16

How to add constraint into loglikelihood function?

1 Answers1