Implementing Neyman Type A distribution in NIMBLE

[Constantin Glen]

Good evening, 

I am trying to implement a model using Neyman's Type A distribution (see:: https://en.m.wikipedia.org/wiki/Neyman_Type_A_distribution). I tried implementing a simple linear model with two parameters (intercept and slope). 

I took the PMF and random number generator from the Wiki page (which seems correct), i.e., 

### load nimble

library(nimble)

### random number generator
rNeymanA <- nimbleFunction(
  run = function(n = integer(0, default = 1),
                 lambda = double(0), phi = double(0)){
    if(n != 1) print("rnorm_vec only allows n = 1; using n = 1.")
    if(lambda < 0  | is.na(lambda)){
      return( NA )
    }else if(phi < 0 | is.na(phi)){
      return( NA )
    }else{
      k <- rpois(1,lambda)
      X <- sum(rpois(k, phi))
    }
    return(X)
    returnType(double(0))
  }
)
compileNimble(rNeymanA)

### PMF for Neyman's type A distribution
dNeymanA <- nimbleFunction(
  run = function(x = double(0), lambda = double(0), phi = double(0),
                 log = integer(0, default = 1)){
    returnType(double(0))
    p    <- nimNumeric(x)
    p[1] <- exp(-lambda + lambda*exp(-phi))
    c    <- lambda*phi*exp(-phi)
    if(x == 0){
      if(log){ return(log(p[1])) }else { return(p[1]) }
    }else{
      for(i in 1:x) {
        sumA = 0
        for (r in 0:(i-1)) {
          sumA = sumA + (phi^r)/( factorial(r) )*p[i-r]
        }
        p[i+1] = (c/(i))*sumA
      }
      Prob <- p[i+1] 
    }
    if(is.na(Prob) | is.nan(Prob)){
      return(-Inf)
    }else{
      if(log){ return(log(Prob)) } else{ return(Prob) }  
    }
  }
)

compileNimble(dNeymanA)
assign('dNeymanA', dNeymanA, envir = .GlobalEnv)
registerDistributions("dNeymanA")

When I try to run the model in nimble, my R crashes and reports a fatal error. However, using maximum likelihood seems to work.

Example: 

### nimble function

NeymanADC <- nimbleCode({
 
  # Prior for delta parameter of Neyman type A distribution
  logphi    ~ dnorm(mean=0, tau=0.0001)
  beta[1:p] ~ dmnorm(zeros[1:p], prec=omega[1:p, 1:p])
  phi       <- exp(logphi)
 
  # likelihood
  for(i in 1:N){
    mu[i] <- exp( beta[1] + beta[2] * X[i] )
    Y[i] ~ dNeymanA(lambda = mu[i]/phi, phi = phi)
  }
})

### using maximum likelihood 

NeymanAnllcov = function(data, par){
  # Data
  Y <- data[,1]
  X <- data[,2]
  # Parameters
  len.par <- length(par)
  mu      <- exp( par[1] + par[len.par-1] * X )
  phi     <- exp( par[len.par] )
  # Likelihood
  llike.vec = sapply(1:length(Y), function(i) dNeymanA( x=Y[i], lambda=mu[i]/phi, phi=phi, 1))
  ret = sum(llike.vec)
  return(-ret)
}

## parameters and data generation
pars <- c(1.5, 0.24, 0.7)
phi  <- exp(pars[3])
nobs <- 500
x1 <- runif(nobs, min=0, max=10)
xb <- pars[1] + pars[2]*x1
exb <- exp(xb)
nApy <- c()
for (i in 1:nobs){
  nApy[i] <- rNeymanA(1, lambda=exb[i]/phi, phi = phi)
}
nAdata <- data.frame(nApy, x1)

### optimize the likelihood function using nlminb

opt <- nlminb(start=c(1,0.1,0.1), objective=NeymanAnllcov,
              control=list(trace=1), data=nAdata)

opt$par

This produces the expected results (estimates close to the true values).

p <- 2; n <- nobs
constants <- list(N = n, K=K, p=p, zeros = rep(0, p), omega = 0.0001 * diag(p))
data <- list(Y = nAdata$nApy, X = nAdata$x1)

tictoc::tic()
example.out <- nimbleMCMC(code = NeymanADC,
                          constants = list(N = n, K=K, p=p, zeros = rep(0, p),
                                           omega = 0.0001 * diag(p)),
                          data = list(Y = nAdata$nApy, X = nAdata$x1),
                          nburnin = 10000, niter = 50000, thin = 5, nchains=5,
                          inits = list(beta = c(5,0.1), phi=0.1),
                          summary = TRUE, samplesAsCodaMCMC=TRUE, WAIC = FALSE,
                          monitors = c("beta", "phi"))
tictoc::toc()

This, however, breaks, and I have no idea where the error is since no error is reported before crashing. I am relatively new to nimble and any advice would help. Thank you. 

Best, 

CGG.

[Perry de Valpine]

Dear Constantin,

This looks interesting and welcome to nimble. I am not 100% sure, but from looking at your code, it appears that in dNeymanA, when i == x at the end of for(i in 1:x), p[i+1] will refer to p[x+1] and that will be beyond the length of p (which was created to have length x). Uncompiled execution is forgiving of this situation because R will automatically extend a vector, but compiled nimble won't do that, and it could result in a crash.

Please let us know if that does not turn out to be the (only) problem.

Also a good strategy for debugging dNeymanA could be to compile it separately,

CdNeymanA <- compileNimble(dNeymanA)

and then call it with various inputs for testing purposes. It is often simpler to figure out a problem this way, without the added layer of the model.

Best wishes and HTH

Perry

--
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 on the web visit https://groups.google.com/d/msgid/nimble-users/a663e405-76d1-4724-8daa-51d03ed7dd95n%40googlegroups.com.