R function produces NAN, Diagnostic. STAT::Optim

Question

I am running this code in order to calculate MLE for Beta Binomial distribution in a Zero-Inflated Model. However I am getting some error for a specific data. Would you please help me out? Here is the R code:

 type = "zi"; lowerbound = 0.01; upperbound = 10000;
    n=27.59554;alpha1=19.22183;alpha2=41.90441;
    x=c(9, 12, 13, 11, 12,  0, 11,  0, 11,  0,  7, 12,  0, 11, 12,  0,  6,  2,  6,  2,  0,  0,  
    9, 10, 10,  0,  0,  0, 10, 11)
    N = length(x)
    t = x[x > 0]  
    m = length(t)
    neg.log.lik <- function(y) 
    { 
      n1 = y[1]
      a1 = y[2]
      b1 = y[3]
      logA = lgamma(a1 + n1 + b1) + lgamma(b1)
      logB = lgamma(a1 + b1) + lgamma(n1 + b1)
      ans = m * log(1 - exp(logB - logA)) + m * logA - m * 
        lgamma(n1 + 1) - sum(lgamma(t + a1)) - sum(lgamma(n1 - 
                  t + b1)) - m * lgamma(a1 + b1) + sum(lgamma(t + 1)) + 
        sum(lgamma(n1 - t + 1)) + m * lgamma(a1)
      return(ans)
    }
    gp <- function(y) 
    {
      #n1=27.59554;a1=19.22183;b1=41.90441;
      n1 = y[1]
      a1 = y[2]
      b1 = y[3]
      logA = lgamma(a1 + n1 + b1) + lgamma(b1)
      logB = lgamma(a1 + b1) + lgamma(n1 + b1)
      dn = -m * exp(logB - logA) * (digamma(n1 + b1) - digamma(a1 + 
                      n1 + b1))/(1 - exp(logB - logA)) - m * digamma(n1 + 
                  1) - sum(digamma(n1 + b1 - t)) + sum(digamma(n1 - 
                                                                                                                                                                 t + 1)) + m * digamma(a1 + n1 + b1)
      da = -m * exp(logB - logA) * (digamma(a1 + b1) - digamma(a1 + 
                  n1 + b1))/(1 - exp(logB - logA)) - sum(digamma(t + 
                   a1)) - m * digamma(a1 + b1) + m * digamma(a1 + n1 + 
                                                                                                                                                              b1) + m * digamma(a1)
      db = -m * exp(logB - logA) * (digamma(a1 + b1) + digamma(n1 + 
                                                                 b1) - digamma(a1 + n1 + b1) - digamma(b1))/(1 - exp(logB - 
                                                                  logA)) + m * digamma(b1) - sum(digamma(n1 - t + b1))-
        m * digamma(a1 + b1) + m * digamma(a1 + n1 + b1)
      return(c(dn, da, db))
    }
    estimate = stats::optim(par = c(n, alpha1, alpha2), fn = neg.log.lik, 
                            gr = gp, method = "L-BFGS-B", lower = c(max(x) - lowerbound, 
                          lowerbound, lowerbound), upper = c(upperbound, upperbound, upperbound))

The error is : Error in stats::optim(par = c(n, alpha1, alpha2), fn = neg.log.lik, gr = gp, : non-finite value supplied by optim In addition: Warning message: In digamma(n1 + b1 - t) : NaNs produced

Do you have any idea? appretiate any suggestions

You need to ensure that digamma gets passed only strictly positive values. The function is not defined for x <= 0. — Roland, Mar 03 '21 at 19:24
I did. t=x[x>0] will do that. I can not figure out what the problem is. — Stat fellow, Mar 03 '21 at 19:53

LocoGris · Answer 1 · 2021-03-03T20:42:15.303

THe issue is that at some point dn becomes NaN due to digamma(0). An option is to take into account that eventuality like I did. But you should explore what to do in that event. You have the issue here digamma(n1 + b1 - t) at some point it produces digamma(0), so NaN so failure.

type = "zi"; lowerbound = 0.01; upperbound = 10000;
n=27.59554;alpha1=19.22183;alpha2=41.90441;
x=c(9, 12, 13, 11, 12,  0, 11,  0, 11,  0,  7, 12,  0, 11, 12,  0,  6,  2,  6,  2,  0,  0,  
    9, 10, 10,  0,  0,  0, 10, 11)
N = length(x)
t = x[x > 0]  
m = length(t)
neg.log.lik <- function(y) 
{ 
  n1 = y[1]
  a1 = y[2]
  b1 = y[3]
  logA = lgamma(a1 + n1 + b1) + lgamma(b1)
  logB = lgamma(a1 + b1) + lgamma(n1 + b1)
  ans = m * log(1 - exp(logB - logA)) + m * logA - m * 
    lgamma(n1 + 1) - sum(lgamma(t + a1)) - sum(lgamma(n1 - 
                                                        t + b1)) - m * lgamma(a1 + b1) + sum(lgamma(t + 1)) + 
    sum(lgamma(n1 - t + 1)) + m * lgamma(a1)
  return(ans)
}
gp <- function(y) 
{
  #n1=27.59554;a1=19.22183;b1=41.90441;
  n1 = y[1]
  a1 = y[2]
  b1 = y[3]
  logA = lgamma(a1 + n1 + b1) + lgamma(b1)

  logB = lgamma(a1 + b1) + lgamma(n1 + b1)

  dn = -m * exp(logB - logA) * (digamma(n1 + b1) - digamma(a1 + 
                                                             n1 + b1))/(1 - exp(logB - logA)) - m * digamma(n1 + 
                                                                                                              1) - sum(digamma(n1 + b1 - t)) + sum(digamma(n1 - 
                                                                                                                                                             t + 1)) + m * digamma(a1 + n1 + b1)
  if(is.na(dn)){
    dn=-99999999
  }
  da = -m * exp(logB - logA) * (digamma(a1 + b1) - digamma(a1 + 
                                                             n1 + b1))/(1 - exp(logB - logA)) - sum(digamma(t + 
                                                                                                              a1)) - m * digamma(a1 + b1) + m * digamma(a1 + n1 + 
                                                                                                                                                          b1) + m * digamma(a1)
  db = -m * exp(logB - logA) * (digamma(a1 + b1) + digamma(n1 + 
                                                             b1) - digamma(a1 + n1 + b1) - digamma(b1))/(1 - exp(logB - 
                                                                                                                   logA)) + m * digamma(b1) - sum(digamma(n1 - t + b1))-
    m * digamma(a1 + b1) + m * digamma(a1 + n1 + b1)
  print("dn")
  print(dn)
  print("da")
  print(da)
  print("db")
  print(db)

  return(c(dn, da, db))
}
estimate = stats::optim(par = c(n, alpha1, alpha2), fn = neg.log.lik, 
                        gr = gp, method = "L-BFGS-B", lower = c(max(x) - lowerbound, 
                                                                lowerbound, lowerbound), upper = c(upperbound, upperbound, upperbound))

Other option is to add a small quantity: digamma(n1 + b1 - t)+0.001

type = "zi"; lowerbound = 0.01; upperbound = 10000;
n=27.59554;alpha1=19.22183;alpha2=41.90441;
x=c(9, 12, 13, 11, 12,  0, 11,  0, 11,  0,  7, 12,  0, 11, 12,  0,  6,  2,  6,  2,  0,  0,  
    9, 10, 10,  0,  0,  0, 10, 11)
N = length(x)
t = x[x > 0]  
m = length(t)
neg.log.lik <- function(y) 
{ 
  n1 = y[1]
  a1 = y[2]
  b1 = y[3]
  logA = lgamma(a1 + n1 + b1) + lgamma(b1)
  logB = lgamma(a1 + b1) + lgamma(n1 + b1)
  ans = m * log(1 - exp(logB - logA)) + m * logA - m * 
    lgamma(n1 + 1) - sum(lgamma(t + a1)) - sum(lgamma(n1 - 
                                                        t + b1)) - m * lgamma(a1 + b1) + sum(lgamma(t + 1)) + 
    sum(lgamma(n1 - t + 1)) + m * lgamma(a1)
  return(ans)
}
gp <- function(y) 
{
  #n1=27.59554;a1=19.22183;b1=41.90441;
  n1 = y[1]
  a1 = y[2]
  b1 = y[3]
  logA = lgamma(a1 + n1 + b1) + lgamma(b1)
  
  logB = lgamma(a1 + b1) + lgamma(n1 + b1)
  
  dn = -m * exp(logB - logA) * (digamma(n1 + b1) - digamma(a1 + 
                                                             n1 + b1))/(1 - exp(logB - logA)) - m * digamma(n1 + 
                                                                                                              1) - sum(digamma(n1 + b1 - t)+0.001) + sum(digamma(n1 - 
                                                                                                                                                             t + 1)) + m * digamma(a1 + n1 + b1)
  da = -m * exp(logB - logA) * (digamma(a1 + b1) - digamma(a1 + 
                                                             n1 + b1))/(1 - exp(logB - logA)) - sum(digamma(t + 
                                                                                                              a1)) - m * digamma(a1 + b1) + m * digamma(a1 + n1 + 
                                                                                                                                                          b1) + m * digamma(a1)
  db = -m * exp(logB - logA) * (digamma(a1 + b1) + digamma(n1 + 
                                                             b1) - digamma(a1 + n1 + b1) - digamma(b1))/(1 - exp(logB - 
                                                                                                                   logA)) + m * digamma(b1) - sum(digamma(n1 - t + b1))-
    m * digamma(a1 + b1) + m * digamma(a1 + n1 + b1)
  print("dn")
  print(dn)
  print("da")
  print(da)
  print("db")
  print(db)
  
  return(c(dn, da, db))
}
estimate = stats::optim(par = c(n, alpha1, alpha2), fn = neg.log.lik, 
                        gr = gp, method = "L-BFGS-B", lower = c(max(x) - lowerbound, 
                                                                lowerbound, lowerbound), upper = c(upperbound, upperbound, upperbound))

Thank you so much for your answer. I appreciate your help. Best. — Stat fellow, Mar 03 '21 at 23:26

R function produces NAN, Diagnostic. STAT::Optim

1 Answers1