PIG distribution using nimble functions

[Laura Pirâmides]

Hi everyone! I have been using nimble for count over dispersed data and tried using the PIG distribution through nimbleFunctions, however, i got some unstable results that I believe are due to my implementation specifically. I will attach my implementation, if anyone has any idea how to improve it, please feel free to reach out. 

# function for the inverse Gaussian (needed for PIG implementation)
dinvgaussfun <- nimbleFunction(
  run = function(x = double(0), lambda = double(0, default = 1),
                 gamma = double(0,default = 1),
                 log = integer(0, default = 0)) {
    returnType(double(0))
    f = sqrt(gamma/(2*pi*x^3)) * exp(-gamma * (x-lambda)^2/(2*x*lambda^2))
    if(log) return(log(f))
    else return(f)
   
  })

rinvgaussfun <- nimbleFunction(
  run = function(n = integer(0, default = 1), lambda = double(0, default = 1),
                 gamma = double(0,default = 1)) {
    returnType(double(0))
    if(n != 1) print("rinvgaussfun só permite n = 1; usando n = 1.")
    y =rnorm(1,0,1)^2
    k = lambda + (lambda^2*y/(2*gamma)) -
      (lambda/(2*gamma) * sqrt(4*lambda*gamma*y + lambda^2*y^2))
    z = runif(1,0,1)
   
    if (z <= lambda / (lambda + k)) return (k)
    else return(lambda^2 / k)
  })


dpigfun <- nimbleFunction(
  run = function(x = double(0), mu = double(0, default = 1),
                 phi = double(0, default = 1),
                 log = integer(0, default = 0)) {
    returnType(double(0))
   
    if (x < 0) return(-Inf)  

   # below I have an approximation for the Bessel function, since it can be unstable
   
    nu <- abs(x-0.5)
    z <- sqrt(phi * (phi + 2 * mu))/nu
    sz <- sqrt(1 + z^2)
    t <- 1/sz
    eta <- sz + log(z/(1 + sz))
   
    t2 <- t^2
    u1.t <- (t * (3 - 5 * t2))/24
    u2.t <- t2 * (81 + t2 * (-462 + t2 * 385))/1152
    u3.t <- t * t2 * (30375 + t2 * (-369603 + t2 * (765765 - t2 * 425425)))/414720
    t4 <- t2 * t2
    u4.t <- t4 * (4465125 + t2 * (-94121676 + t2 * (349922430 + t2 * (-446185740 + t2 * 185910725))))/39813120
    u5.t <- t * t4 * (1519035525 + t2 * (-49286948607 + t2 * (284499769554 + t2 * (-614135872350 + t2 * (566098157625 - t2 * 188699385875)))))/6688604160
    d <- (-u1.t + (u2.t + (-u3.t + (u4.t - u5.t/nu)/nu)/nu)/nu)/nu
   
    bessel_log <- log1p(d) - nu * eta - (log(sz) - log(pi/(2 * nu)))/2
   
    #   density
    log_f <- 0.5 * log(2 / pi) + phi + x * log(mu * phi) - lgamma(x + 1) -
      0.5 * (x - 0.5) * log(phi * (phi + 2 * mu)) + bessel_log
   
   
     if (log) return(log_f)
    else return(exp(log_f))
  })

rpigfun <- nimbleFunction(
  run = function(n = double(0, default = 1), mu = double(0, default = 1),
                 phi = double(0, default = 1)) {
    returnType(double(0))
   
    if (n != 1) print("rpigfun só permite n = 1; usando n = 1.")
   
    # z ~ Inverse Gaussian
    z <- rinvgaussfun(1, lambda = 1, gamma = phi)  
   
    return(rpois(1, z * mu))
  })

Also, I have implemented it using the Bessel function available in nimble, this approximation was done as a way to minimize the instability.

[Chris Paciorek]

Hi Laura,

In general, it's best to calculate density values on the log scale (and in fact, nimble works with densities directly on the log scale).

So your dinvgaussfun would look like this (check my math!):

    logf <- 0.5*log(lambda) - log(pi) -3*log(x) - gamma * (x-lambda)^2/(2*x*lambda^2))

    if(log) return(logf)
    else return(exp(f))

I'm not sure exactly what you mean by unstable, so I don't know if the above change will fix the problem. If not, please provide more details/an example of the instability.

And regardless of whether the above change fixes the problem, I am curious what instability you see with our besselK function. That uses a function from the R math library, which I would think would be robustly implemented.

-chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/09c25984-9b74-4596-9e7f-db7bdcfda5e7n%40googlegroups.com.