Nimble MCMC algorithm no working with Conway Maxwell Poisson Distribution
122 views
Skip to first unread message
Lokesh Arya
Mar 31, 2019, 5:24:33 PM3/31/19
to nimble-users
Hi All,
I'm trying to estimate the parameters of Zero-Inflated Conway Maxwell Poisson distribution in Nimble using MCMC algorithm.
I took reference from the example provided on the Nimble website for Zero-inflated Poisson distribution.
I was expecting the same behavior of the algorithm as in Zero-Inflated Poisson example, But I don't why it is behaving so weirdly,
Can somebody please help me out in this situation.
library(CompGLM)
library(nimble)
ZICMPcode <- nimbleCode({
p~dunif(0,1)
lam~dunif(0,100)
nu~dunif(0,100)
for(i in 1:N) {
y[i] ~ dZICMP(lam, nu,sumTo = 200L,zeroProb = p)
}
} )
## constants, data, and initial values
constants <- list(N = 100)
dCMP <- nimbleRcall(function(y = integer(), lam = double(), nu = double(),sumTo = integer(),logP = logical(0, default = 0)){}, Rfun = 'dcomp',
returnType = double(), envir = .GlobalEnv)
rCMP <- nimbleRcall(function(n = integer(), lam = double(), nu = double(),sumTo = integer()){}, Rfun = 'rcomp',
returnType = integer(), envir = .GlobalEnv)
dZICMP <- nimbleFunction(
run = function(x = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double(), log = logical(0, default = 0)) {
returnType(double())
## First handle non-zero data
if(x != 0) {
## return the log probability if log = TRUE
if(log) return(dCMP(x, lam, nu,sumTo = 200L, log = TRUE) + log(1-zeroProb))
## or the probability if log = FALSE
else return((1-zeroProb) * dCMP(x, lam, nu,sumTo = 200L, log = FALSE))
}
## From here down we know x is 0
totalProbZero <- zeroProb + (1-zeroProb) * dCMP(0, lam, nu,sumTo = 200L, log = FALSE)
if(log) return(log(totalProbZero))
return(totalProbZero)
})
rZICMP <- nimbleFunction(
run = function(n = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double()) {
returnType(integer())
isStructuralZero <- rbinom(1, prob = zeroProb, size = 1)
if(isStructuralZero) return(0)
return(rCMP(1, lam,nu,sumTo = 200L))
})
registerDistributions(list(
dZICMP = list(
BUGSdist = "dZICMP(lam, nu,sumTo, zeroProb)",
discrete = TRUE,
range = c(0, Inf),
types = c('value = integer()', 'lam = double()', 'nu = double()','sumTo = integer()', 'zeroProb = double()')
)))
ZICMPmodel <- nimbleModel( ZICMPcode,constants=constants,check = FALSE)
ZICMPmodel$p <- .4
ZICMPmodel$lam <- 1.8
ZICMPmodel$nu<-0.9
ZICMPmodel$simulate('y')
simulatedData <- ZICMPmodel$y
simulatedData
hist(simulatedData)
ZICMPmodel$setData(list(y = simulatedData))
cZICMPmodel <- compileNimble(ZICMPmodel)
ZICMPmcmc <- buildMCMC(ZICMPmodel)
cZICMPmcmc <- compileNimble(ZICMPmcmc, project = ZICMPmodel)
cZICMPmcmc$run(10000)
samples <- as.matrix(cZICMPmcmc$mvSamples)
summary(samples)
plot(samples[,'lam'], main = 'lambda trace plot')
plot(samples[,'nu'], main = 'nu trace plot')
plot(samples[,'p'], main = 'zeroProb trace plot')
Thank you
Daniel Turek
Apr 1, 2019, 10:15:47 AM4/1/19
to Lokesh Arya, nimble-users
Lokesh, thanks for your use of NIMBLE. I'm not quite sure what the problem you were having was. I provided initial values to the model, so that the first model build was "fully initialized", which I think is always a good idea. Then I continued with your mechanism of generating data (which turns out to be the same: all 1's), and using that.
Then, I also broke down the MCMC building into the separate steps of configuration, then building, so we can inspect the MCMC itself. I think this is *always* a useful intermediate step, so you know what's actually going on in the MCMC, and can catch if something (usually in the samplers) seems amiss.
Everything compiled and ran fine. You'll definitely want to do longer MCMC runs (I lowered it to 1000 iterations), but the MCMC seems to get started fine.
You might consider using the runMCMC() function, rather than mcmc$run(), but either will work fine.
You might also check out the samplesPlot() function, provided in the package basicMCMCplots. I used it in the code below.
Let us know if you're still having problems here. All my code is copied below. Most of my comments that I added are proceeded by "DT"
Cheers!
Daniel
library(CompGLM)
library(nimble)
ZICMPcode <- nimbleCode({
p~dunif(0,1)
lam~dunif(0,100)
nu~dunif(0,100)
for(i in 1:N) {
y[i] ~ dZICMP(lam, nu,sumTo = 200L,zeroProb = p)
}
})
## constants, data, and initial values
constants <- list(N = 100)
dCMP <- nimbleRcall(function(y = integer(), lam = double(), nu = double(),sumTo = integer(),logP = logical(0, default = 0)){}, Rfun = 'dcomp',
returnType = double() )## DT: shouldn't need this: envir = .GlobalEnv)
rCMP <- nimbleRcall(function(n = integer(), lam = double(), nu = double(),sumTo = integer()){}, Rfun = 'rcomp',
returnType = integer()) ## DT: shouldn't need this: envir = .GlobalEnv)
## DT: let's include the model check.
## it's helpful, especially for such a small model:
ZICMPmodel <- nimbleModel( ZICMPcode, constants=constants) ##check = FALSE)
ZICMPmodel$calculate() ## something is wrong here!
## DT:
## Error in if (x != 0) { : missing value where TRUE/FALSE needed
## looks like since we're missing values for y[i] (no initial values were provided)
## this first logical if() statement in your dZICMP() function definition
## is failing. Let's provide initial values and try again:
inits <- list(p = 0.5,
lam = 1,
nu = 1)
data <- list(y = rep(1, constants$N))
ZICMPmodel <- nimbleModel( ZICMPcode, constants=constants, data=data, inits=inits)
ZICMPmodel$calculate() ## -178.5251 DT: that's better
ZICMPmodel$p <- .4
ZICMPmodel$lam <- 1.8
ZICMPmodel$nu<-0.9
ZICMPmodel$simulate('y')
ZICMPmodel$calculate() ## -191.4907 DT: still good
simulatedData <- ZICMPmodel$y
simulatedData
hist(simulatedData)
ZICMPmodel$setData(list(y = simulatedData))
cZICMPmodel <- compileNimble(ZICMPmodel)
## DT: let's configure the MCMC first, look at the samplers
## that are assigned, and the monitors. Always
## helpful and a good idea to know what's going on in the MCMC:
##ZICMPmcmc <- buildMCMC(ZICMPmodel) not this line, but instead:
conf <- configureMCMC(ZICMPmodel)
conf$printMonitors()
## thin = 1: p, lam, nu
conf$printSamplers()
## [1] RW sampler: p
## [2] RW sampler: lam
## [3] RW sampler: nu
## DT: now build the MCMC, and compile it:
ZICMPmcmc <- buildMCMC(conf)
cZICMPmcmc <- compileNimble(ZICMPmcmc, project = ZICMPmodel)
cZICMPmcmc$run(1000) ## using fewer iterations
samples <- as.matrix(cZICMPmcmc$mvSamples)
samplesSummary(samples)
## Mean Median St.Dev. 95%CI_low 95%CI_upp
## lam 0.62186699 0.560347129 0.18223126 0.4068033525 0.93935815
## nu 0.34242280 0.014943882 0.40012792 0.0026253462 0.90000000
## p 0.02232572 0.002589668 0.02716132 0.0009254903 0.05182992
library(basicMCMCplots)
samplesPlot(samples) ## DT: looks like we need longer runs, but things are off to a good start
--
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 post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/6a5fdefb-261c-4323-b8d0-2426b5f05f41%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Lokesh Arya
Apr 12, 2019, 7:37:21 AM4/12/19
to nimble-users
Thanks for the help Daniel.
but it is still not working.
actually, my main objective is to make a Zero-inflated Conway Maxwell Poisson model through MCMC
And I'm don't know why I'm getting a pattern in the sampleplot for the coefficients.
can you please help me out.
Thank you again
library(CompGLM)
library(nimble)
x1<-rnorm(100)
lambda<-exp(0.5+0.9*x1)
nu<-exp(1+0.8*x1)
p<-exp(0.3+0.7*x1)/(1+exp(0.3+0.7*x1))
bin<-c()
for(i in 1:100)
bin[i] = rbinom(1,1,p[i])
obs<-c()
for(i in 1:100){
obs[i] = ifelse(bin[i] == 0,0,rcomp(1,lambda[i],nu[i]))
}
hist(obs)
ZICMPcode <- nimbleCode({
alpha0 ~ dnorm(0,sd = 10000)
alpha1 ~ dnorm(0,sd = 10000)
beta0 ~ dnorm(0,sd = 10000)
beta1 ~ dnorm(0,sd = 10000)
gamma0 ~ dnorm(0,sd = 10000)
for(i in 1:N) {
logit(p[i])<-alpha0+alpha1*x[i]
log(lam[i])<-beta0+beta1*x[i]
log(nu[i])<-gamma0
y[i] ~ dZICMP(lam[i], nu[i],sumTo = 200L,zeroProb = p[i])## Note NIMBLE allows R-like named-parameter syntax
}
} )
## constants, data, and initial values
constants <- list(N = 100)
data <- list(y = obs,x = x1)
inits <- list(alpha0 = 0,alpha1 = 0,beta0 = 0, beta1 = 0,gamma0=0)
ZICMPmodel <- nimbleModel( ZICMPcode,constants=constants,data = data,inits = inits,check = T)
ZICMPmodel$calculate()
#ZICMPmodel$alpha0 <- .4 ## Choose values of p and lambda
#ZICMPmodel$alpha1 <- 1.8
#ZICMPmodel$beta0<-0.5
#ZICMPmodel$beta1<-0.5
#ZICMPmodel$gamma0<-0.5
#ZICMPmodel$simulate('y') ## Simulate values of y[1]...y[100]
#simulatedData <- ZICMPmodel$y
#simulatedData
#hist(simulatedData)
#ZICMPmodel$setData(list(y = simulatedData)) ## Set those values as data in the model
ZICMPmodel$calculate()
cZICMPmodel <- compileNimble(ZICMPmodel) ## Compile the model
conf <- configureMCMC(ZICMPmodel)
conf$printMonitors()
conf$printSamplers()
ZICMPmcmc <- buildMCMC(conf)
cZICMPmcmc <- compileNimble(ZICMPmcmc, project = ZICMPmodel)
cZICMPmcmc$run(10000)
samples <- as.matrix(cZICMPmcmc$mvSamples)
samplesSummary(samples)
library(basicMCMCplots)
samplesPlot(samples)
To unsubscribe from this group and stop receiving emails from it, send an email to nimble...@googlegroups.com.
Lokesh Arya
Apr 12, 2019, 7:40:14 AM4/12/19
to nimble-users
please note the correction
nu<- exp(1)
Perry de Valpine
Apr 12, 2019, 11:01:00 AM4/12/19
to Lokesh Arya, nimble-users
When I replace
nu<-exp(1+0.8*x1)
with
nu<- exp(1)
as you stated, the simulated obs are all 0 or nearly all 0. Is that part of the problem?
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/c9fdde91-1074-4d4f-ae6b-8646cf5c94dd%40googlegroups.com.
Lokesh Arya
Apr 12, 2019, 11:08:37 AM4/12/19
to nimble-users
No that's not the problem
I'm intentionally doing that to simulate Zero-Inflated Data.
problem is that
"The parameters are not converging, they are showing an upward and downward trend kind of pattern instead of an expected random pattern as shown in the Zero-inflated Poisson example on the website."
Lokesh Arya
Apr 12, 2019, 12:28:34 PM4/12/19
to nimble-users
I'm trying to simulate the following kind of situation and wanna fit a Zero-inflated Conway Maxwell Poisson model over it using MCMC algorithm
set.seed(124)
x1<-rnorm(100)
lambda<-exp(3.5+0.9*x1)
nu<-exp(1)
p<-plogis(0.3+0.7*x1)
bin<-c()
for(i in 1:100)
bin[i] = rbinom(1,1,p[i])
obs<-c()
for(i in 1:100){
obs[i] = ifelse(bin[i] == 0,0,rcomp(1 , lam = lambda[i],nu = nu))
}
hist(obs)
i'm getting output(sample plot file) after running MCMC algorithm
Perry de Valpine
Apr 12, 2019, 2:19:30 PM4/12/19
to Lokesh Arya, nimble-users
Hi Lokesh,
Thanks for making an example with more informative data, so it's more clear that there is a real sampling problem.
I think this is revealing a bug related to nimbleRcall and R's garbage collection, which is creating incorrect values being passed from C++ to R or back. We have seen a bug related to this previously, which we fixed, but evidently there is more to it. I think it has to do with use of PROTECT in C++ code for R.
What convinces me that something like this is the problem is that I was able to write dcomp purely as a nimbleFunction, so that no nimbleRcall is needed, and now the results look reasonable. To write it this way, I looked at the CompGLM C++ source code, where I saw that the functions dcomp and Z were both easy to write as nimbleFunctions.
Using the following pieces makes it work for me. I hope this gets you under way on your goal, while we take longer to figure out what the problem with nimbleRcall is.
calcZ <- nimbleFunction(
run = function(lam = double(), nu = double(), sumTo = double()) {
## CompGLM C++ source code:
## double Z(double lam, double nu, int sumTo) {
## double sum = 1.0;
## double factorial = 1.0;
## double lamPower = lam;
## for (int i = 1; i <= sumTo; ++i) {
## factorial *= i;
## sum += lamPower / pow(factorial, nu);
## lamPower *= lam;
## }
## return sum;
## }
sum <- 1
factorial <- 1
lamPower <- lam
for(i in 1:sumTo) {
factorial <- factorial * i
sum <- sum + lamPower / pow(factorial, nu)
lamPower <- lamPower * lam
}
return(sum)
returnType(double())
}
)
dCMP_nf <- nimbleFunction(
## Normally a dist random variable must be called x, but this is not really the dist function being called, so the name "y" will work
run = function(y = double(), lam = double(), nu = double(), sumTo = double(), log = integer(0, default = 0) ) {
## C source code from package CompGLM (which is wrapped by another layer that handles log):
## double dcomp(int y, double lam, double nu, int sumTo) {
## if (y < 0) {
## return 0.0;
## } else {
## return pow(lam, y) / (pow(factorial(y), nu) * Z(lam, nu, sumTo));
## }
## }
returnType(double())
if(y < 0) {
if(log) return(-Inf)
else return(0)
} else {
prob <- pow(lam, y) / (pow(factorial(y), nu) * calcZ(lam, nu, sumTo))
if(log) return(log(prob))
else return(prob)
}
}
)
dZICMP <- nimbleFunction(
run = function(x = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double(), log = logical(0, default = 0)) {
returnType(double())
## print("x = ", x, " lam = ", lam, " nu = ", nu, " zeroProb = ", zeroProb)
## First handle non-zero data
if(x != 0) {
## return the log probability if log = TRUE
if(log) return(dCMP_nf(x, lam, nu,sumTo = 200L, log = TRUE) + log(1-zeroProb))
## or the probability if log = FALSE
else return((1-zeroProb) * dCMP_nf(x, lam, nu,sumTo = 200L, log = FALSE))
}
## From here down we know x is 0
totalProbZero <- zeroProb + (1-zeroProb) * dCMP_nf(0, lam, nu,sumTo = 200L, log = FALSE)
if(log) return(log(totalProbZero))
return(totalProbZero)
})
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/4fd30d5a-1376-452d-a3d5-9579ef5b3ed8%40googlegroups.com.
Lokesh Arya
Apr 12, 2019, 2:48:26 PM4/12/19
to nimble-users
Thank you so much pdevalpine.
It finally worked.
Lokesh Arya
Apr 12, 2019, 6:09:33 PM4/12/19
to nimble-users
Hi Perry
I really appreciate your time and help.
But just wanna ask you one last question.
I'm further trying to apply MCEM algorithm on the above code with adding a random effect in lambda.
mles<-ZICMPMCEM$run(initM = 10000)
warning: logProb of data node y[62]: logProb is NA or NaN.
warning: logProb of data node y[89]: logProb is NA or NaN.
warning: logProb of data node y[62]: logProb is NA or NaN.
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
Error in optim(par = theta, fn = cCalc_E_llk$run, oldParamValues = thetaPrev, :
initial value in 'vmmin' is not finite
I don't why it is showing these warnings and in the end, it is giving me error like this.
Also, it would be very great of you if can check why rCMP_nf function not working, I took that directly from compGLM library C++ code.
right now I put the rCMP_nf function in comment.
library(CompGLM)
library(nimble)
set.seed(124)
x1<-rnorm(100)
x2<-rnorm(100)
lambda<-exp(3.5+0.9*x1+x2)
nu<-exp(1)
p<-plogis(0.3+0.7*x1)
bin<-c()
for(i in 1:100)
bin[i] = rbinom(1,1,p[i])
obs<-c()
for(i in 1:100){
obs[i] = ifelse(bin[i] == 0,0,rcomp(1 , lam = lambda[i],nu = nu))
}
hist(obs)
ZICMPcode <- nimbleCode({
alpha0 ~ dnorm(0,sd = 10000)
alpha1 ~ dnorm(0,sd = 10000)
beta0 ~ dnorm(0,sd = 10000)
beta1 ~ dnorm(0,sd = 10000)
gamma0 ~ dnorm(0,sd = 10000)
sigma~dunif(0,100)
for(i in 1:N) {
u[i]~dnorm(0,sigma)
logit(p[i])<-alpha0+alpha1*x[i]
log(lam[i])<-beta0+beta1*x[i]+u[i]
log(nu[i])<-gamma0
y[i] ~ dZICMP(lam[i], nu[i],sumTo = 200L,zeroProb = p[i])## Note NIMBLE allows R-like named-parameter syntax
}
} )
## constants, data, and initial values
constants <- list(N = 100)
data <- list(y = obs,x = x1)
inits <- list(alpha0 = rnorm(1),alpha1 = rnorm(1),beta0 = rnorm(1), beta1 = rnorm(1),gamma0=rnorm(1))
#dCMP <- nimbleRcall(function(y = integer(), lam = double(), nu = double(),sumTo = integer(),logP = logical(0, default = 0)){}, Rfun = 'dcomp',
# returnType = double(), envir = .GlobalEnv)
rCMP <- nimbleRcall(function(n = integer(), lam = double(), nu = double(),sumTo = integer()){}, Rfun = 'rcomp',
returnType = integer(), envir = .GlobalEnv)
#rCMP_nf<-nimbleFunction(
# run = function(n = integer(),lam = double(),nu = double(),sumTo = double()){
# returnType(integer())
# dist<-c()
# for(i in 1:sumTo){
# dist[i]<-dCMP_nf(i,lam,nu,sumTo)
# }
# sample<-c()
# for(i in 1:n){
# u=runif(0,1)
# sum = 0
# for(j in 1:sumTo){
# sum = sum + dist[j]
# if(u<sum){
# sample[i] = j
# break
# }
# }
# }
# return(sample)
# }
#)
dZICMP <- nimbleFunction(
run = function(x = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double(), log = logical(0, default = 0)) {
returnType(double())
## print("x = ", x, " lam = ", lam, " nu = ", nu, " zeroProb = ", zeroProb)
## First handle non-zero data
if(x != 0) {
## return the log probability if log = TRUE
if(log) return(dCMP_nf(x, lam, nu,sumTo = 200L, log = TRUE) + log(1-zeroProb))
## or the probability if log = FALSE
else return((1-zeroProb) * dCMP_nf(x, lam, nu,sumTo = 200L, log = FALSE))
}
## From here down we know x is 0
totalProbZero <- zeroProb + (1-zeroProb) * dCMP_nf(0, lam, nu,sumTo = 200L, log = FALSE)
if(log) return(log(totalProbZero))
return(totalProbZero)
})
rZICMP <- nimbleFunction(
run = function(n = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double()) {
returnType(integer())
isStructuralZero <- rbinom(1, prob = zeroProb, size = 1)
if(isStructuralZero) return(0)
return(rCMP(1, lam,nu,sumTo = 200L))
})
registerDistributions(list(
dZICMP = list(
BUGSdist = "dZICMP(lam, nu,sumTo, zeroProb)",
discrete = TRUE,
range = c(0, Inf),
types = c('value = integer()', 'lam = double()', 'nu = double()','sumTo = integer()', 'zeroProb = double()')
)))
ZICMPmodel <- nimbleModel( ZICMPcode,constants=constants,data = data,inits = inits,check = F)
ZICMPmodel$calculate()
cZICMPmodel <- compileNimble(ZICMPmodel) ## Compile the model
conf <- configureMCMC(ZICMPmodel,nodes = NULL)
conf$printMonitors()
conf$addSampler(type = 'RW', target ='alpha0')
conf$addSampler(type = 'RW', target ='alpha1')
conf$addSampler(type = 'RW', target ='beta0')
conf$addSampler(type = 'RW', target ='beta1')
conf$addSampler(type = 'RW', target ='gamma0')
conf$addSampler(type = 'RW', target ='sigma')
conf$printSamplers()
#ZICMPmcmc <- buildMCMC(conf)
#cZICMPmcmc <- compileNimble(ZICMPmcmc, project = ZICMPmodel)
#sample<-runMCMC(cZICMPmcmc,thin = 1, niter = 10000, nburnin = 2000, nchains = 1, setSeed = TRUE)
#samplesSummary(sample)
#library(basicMCMCplots)
#samplesPlot(sample)
ZICMPmodel2 <- ZICMPmodel$newModel()
box = list( list("alpha0",c(-Inf, Inf),
list("alpha1",c(-Inf, Inf)),
list("beta0",c(-Inf, Inf)),
list("beta1",c(-Inf, Inf)),
list("gamma0"),c(-Inf, Inf)),
list("sigma",c(0,100)))
ZICMPMCEM <- buildMCEM(model = ZICMPmodel2,latentNodes = ZICMPmodel2$getNodeNames(latentOnly = TRUE),burnIn = 1000,mcmcControl = list(adaptInterval = 10), boxConstraints = box,buffer = 1e-6)
mles<-ZICMPMCEM$run(initM = 10000)
Thank you so much again
Perry de Valpine
Apr 12, 2019, 8:35:44 PM4/12/19
to Lokesh Arya, nimble-users
I don't think you need rCMP_nf to work. If you want it to work, a quick guess is that "break" is the problem. When you define the nimbleFunction, there is a warning message about break. That is telling you that it may not compile unless it is going to be defined as a nimbleFunction later (and actually, in this case, it is not even a function call). I think you could replace "break" with "return(sample)".
For MCEM, you may need to provide box constraints for the optimization(s). You could run it uncompiled and trap the warnings or inspect the model variables and log probabilities after it breaks. Or the model might just not be initialized well. Did you try calculating the whole model? You can access model variables and log probabilities and calculate parts of the model in order to track down where there is a problem.
It's a hasty answer at the end of the day in my time zone, but I hope it helps.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/d1d395fd-08b3-46de-af93-45aa47f4f37d%40googlegroups.com.
Lokesh Arya
Apr 14, 2019, 6:52:40 AM4/14/19
to nimble-users
Thank for the help Perry.
MCEM earlier problem resolved but now it showing optimization error
Error in optim(par = theta, fn = cCalc_E_llk$run, oldParamValues = thetaPrev, :
non-finite finite-difference value [4]
Googled it, Peoples are suggesting try different initial values and i did that too but i'm still getting this error. In fact, algo runs for 5-6 iterations also and after that it gives error like this.
there any way to resolve this problem by any way?
Thank you again.
calZ<-function(lam, nu, sumTo){
sum <- 1
factorial <- 1
lamPower <- lam
for(i in 1:sumTo) {
factorial <- factorial * i
sum <- sum + pow((lamPower /factorial), nu)
lamPower <- lamPower * lam
}
return(sum)
}
dCOMPois<-function(y, lam, nu, sumTo, log){
if(y < 0) {
if(log) return(-Inf)
else return(0)
} else {
prob <- pow((pow(lam, y)/factorial(y)),nu)/calZ(lam, nu, sumTo)
if(log) return(log(prob))
else return(prob)
}
}
rCOMPois<-function(n ,lam ,nu ,sumTo){
dist<-c()
for(i in 1:sumTo){
dist[i]<-dCOMPois(i,lam,nu,log =F,sumTo)
}
sample<-c()
for(i in 1:n){
u<-runif(1,0,1)
sum = 0
for(j in 1:sumTo){
sum = sum + dist[j]
if(u < sum){
sample[i] = j
break
}
}
}
return(sample)
}
library(nimble)
set.seed(124)
x1<-rnorm(100)
x2<-rnorm(100)
lambda<-exp(3.5+0.9*x1+x2)
nu<-0.4
p<-plogis(0.3+0.7*x1)
bin<-c()
for(i in 1:100)
bin[i] = rbinom(1,1,p[i])
obs<-c()
for(i in 1:100){
obs[i] = ifelse(bin[i] == 1,0,rCOMPois(1 , lam = lambda[i],nu = nu ,sumTo = 101 ))
}
hist(obs)
ZICMPcode <- nimbleCode({
alpha0 ~ dnorm(0,sd = 100);
alpha1 ~ dnorm(0,sd = 100);
beta0 ~ dnorm(0,sd = 100);
beta1 ~ dnorm(0,sd = 100);
sigma ~ dunif(00,100)
nu ~ dunif(0,100);
for(i in 1:N) {
u[i] ~ dnorm(0,sigma)
logit(p[i])<-alpha0+alpha1*x[i]
log(lam[i])<-beta0+beta1*x[i]+u[i]
y[i] ~ dZICMP(lam[i], nu,sumTo = 101,zeroProb = p[i])## Note NIMBLE allows R-like named-parameter syntax
}
}
)
## constants, data, and initial values
constants <- list(N = 100,x = x1)
data <- list(y = obs)
inits <- list(alpha0 = rnorm(1),alpha1 = rnorm(1),beta0 = rnorm(1), beta1 = rnorm(1) ,nu = rnorm(1),sigma = 1,u = rep(0,100) )
rCMP <- nimbleRcall(function(n = integer(), lam = double(), nu = double(),sumTo = integer()){}, Rfun = 'rCOMPois',
returnType = integer(), envir = .GlobalEnv)
calcZ <- nimbleFunction(
run = function(lam = double(), nu = double(), sumTo = double()) {
## CompGLM C++ source code:
## double Z(double lam, double nu, int sumTo) {
## double sum = 1.0;
## double factorial = 1.0;
## double lamPower = lam;
## for (int i = 1; i <= sumTo; ++i) {
## factorial *= i;
## sum += lamPower / pow(factorial, nu);
## lamPower *= lam;
## }
## return sum;
## }
sum <- 1
factorial <- 1
lamPower <- lam
for(i in 1:sumTo) {
factorial <- factorial * i
sum <- sum + pow((lamPower /factorial), nu)
lamPower <- lamPower * lam
}
return(sum)
returnType(double())
}
)
dCMP_nf <- nimbleFunction(
## Normally a dist random variable must be called x, but this is not really the dist function being called, so the name "y" will work
run = function(y = double(), lam = double(), nu = double(), sumTo = double(), log = integer(0, default = 0) ) {
## C source code from package CompGLM (which is wrapped by another layer that handles log):
## double dcomp(int y, double lam, double nu, int sumTo) {
## if (y < 0) {
## return 0.0;
## } else {
## return pow(lam, y) / (pow(factorial(y), nu) * Z(lam, nu, sumTo));
## }
## }
returnType(double())
if(y < 0) {
if(log) return(-Inf)
else return(0)
} else {
prob <- pow((pow(lam, y)/factorial(y)),nu)/calcZ(lam, nu, sumTo)
if(log) return(log(prob))
else return(prob)
}
}
)
dZICMP <- nimbleFunction(
run = function(x = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double(), log = logical(0, default = 0)) {
returnType(double())
## print("x = ", x, " lam = ", lam, " nu = ", nu, " zeroProb = ", zeroProb)
## First handle non-zero data
if(x != 0) {
## return the log probability if log = TRUE
if(log) return(dCMP_nf(x, lam, nu,sumTo = 101, log = TRUE) + log(1-zeroProb))
## or the probability if log = FALSE
else return((1-zeroProb) * dCMP_nf(x, lam, nu,sumTo = 101, log = FALSE))
}
## From here down we know x is 0
totalProbZero <- zeroProb + (1-zeroProb) * dCMP_nf(0, lam, nu,sumTo = 101, log = FALSE)
if(log) return(log(totalProbZero))
return(totalProbZero)
})
rZICMP <- nimbleFunction(
run = function(n = integer(), lam = double(), nu = double(),sumTo = integer(), zeroProb = double()) {
returnType(integer())
isStructuralZero <- rbinom(1, prob = zeroProb, size = 1)
if(isStructuralZero) return(0)
return(rCMP(1, lam,nu,sumTo = 101))
})
registerDistributions(list(
dZICMP = list(
BUGSdist = "dZICMP(lam, nu,sumTo, zeroProb)",
discrete = TRUE,
range = c(0, Inf),
types = c('value = integer()', 'lam = double()', 'nu = double()','sumTo = integer()', 'zeroProb = double()')
)))
ZICMPmodel <- nimbleModel( ZICMPcode,constants=constants,data = data,inits = inits,check = F)
ZICMPmodel$calculate()
box = list( list("alpha0",c(0, 2)),
list("alpha1",c(0, 2)),
list("beta0",c(0, 5)),
list("beta1",c(0, 2)),
list("nu",c(0, 1)),
list("sigma",c(0,Inf)))
ZICMPMCEM <- buildMCEM(model = ZICMPmodel,latentNodes ='u[1:100]',burnIn = 100, boxConstraints = box,buffer = 1e-6)
mles<-ZICMPMCEM$run(initM = 1000)
Chris Paciorek
Apr 15, 2019, 10:11:27 PM4/15/19
to Lokesh Arya, nimble-users
Hi Lokesh,
This same error message came up with another user in a recent post to
nimble-users. In that case the issue was
that for some set of values of the latent parameters and parameters
being optimized over, the log-likelihood calculation
returns NA or NaN or Inf or -Inf. At this point the expected
log-likelihood is not a finite numeric value and the optimization step
fails.
So my high-level suggestion is that you look and see if there are
latent parameter values that might be simulated
in the embedded MCMC or hyperparameter values that might be used by
the optimizer that would cause an NA/NaN/Inf/-Inf.
In particular think about parameter values on the boundary of finite
parameter ranges. It's possible that because of
numerical round-off error, the MCEM could produce parameter values
that to the computer are exactly on the boundary.
Another thing to look at is that in both calcZ and dCMP_nf you should
probably work on the log scale as much
as possible for better numerical performance - multiplying and
dividing numbers can lead to numerical
overflow and underflow.
E.g. instead of:
prob <- pow((pow(lam, y)/factorial(y)),nu)/calcZ(lam, nu, sumTo)
if(log) return(log(prob))
else return(prob)
calculate the log prob directly and only exponentiate if the "log"
argument is TRUE.
Ideally NIMBLE's MCEM would do a better job trapping this situation
and returning helpful error message...
> To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
This same error message came up with another user in a recent post to
nimble-users. In that case the issue was
that for some set of values of the latent parameters and parameters
being optimized over, the log-likelihood calculation
returns NA or NaN or Inf or -Inf. At this point the expected
log-likelihood is not a finite numeric value and the optimization step
fails.
So my high-level suggestion is that you look and see if there are
latent parameter values that might be simulated
in the embedded MCMC or hyperparameter values that might be used by
the optimizer that would cause an NA/NaN/Inf/-Inf.
In particular think about parameter values on the boundary of finite
parameter ranges. It's possible that because of
numerical round-off error, the MCEM could produce parameter values
that to the computer are exactly on the boundary.
Another thing to look at is that in both calcZ and dCMP_nf you should
probably work on the log scale as much
as possible for better numerical performance - multiplying and
dividing numbers can lead to numerical
overflow and underflow.
E.g. instead of:
prob <- pow((pow(lam, y)/factorial(y)),nu)/calcZ(lam, nu, sumTo)
if(log) return(log(prob))
else return(prob)
argument is TRUE.
Ideally NIMBLE's MCEM would do a better job trapping this situation
and returning helpful error message...
> To post to this group, send email to nimble...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/c535bf0c-a999-43b4-a1d4-4826cc276d87%40googlegroups.com.
Lokesh Arya
Apr 15, 2019, 11:29:02 PM4/15/19
to Chris Paciorek, nimble-users
Thank you so much chris
the problem of getting -Inf is resolved by using log scale
but problem ofoptimization is still happening.
can we do any thing about this?
Thanks again for your help
Chris Paciorek
Apr 16, 2019, 10:22:01 AM4/16/19
to Lokesh Arya, nimble-users
I'm not sure what you mean by "problem of optimization is still
happening". Do you mean that
you are still getting the error "non-finite finite-difference value"
from the optimization?
Did you determine that there are not parameter settings that would
lead to Inf/-Inf/NA/NaN
values in the log-likelihood?
I might be able to take a look in more detail in the next week, but I
can't promise I will have time.
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAJhvOVMDEa1sWMVFCQzUmRDnZzCYaMNZ-fL2nyg%2BP7wLEGVWkA%40mail.gmail.com.
happening". Do you mean that
you are still getting the error "non-finite finite-difference value"
from the optimization?
Did you determine that there are not parameter settings that would
lead to Inf/-Inf/NA/NaN
values in the log-likelihood?
I might be able to take a look in more detail in the next week, but I
can't promise I will have time.
Chris Paciorek
Apr 29, 2019, 10:51:51 AM4/29/19
to Lokesh Arya, nimble-users
Hi Lokesh,
I was finally able to get a bit of time to look more at your MCEM issue. It is indeed a numerical problem caused by not working on the log scale.
It's occurring when the optimization step puts values from the MCMC for the latent nodes into the model for the current set of values for the top level nodes
and the resulting model log probability is Inf or -Inf or NaN. In your case, because calcZ was not being done on the log scale, it was returning Inf or -Inf.
I've attached a version of your code (the code from your April 14 email) where key steps are done on the log scale (see everything where I made a comment noting "CHANGED").
The MCEM for the revised code now runs without the error. However, it takes a very long time to run/converge.
I'm also going to discuss with our development team some ideas so that when this sort of thing happens, we catch the Inf/-Inf/NaN that occurs in the calculation and
give the user more information about what has gone wrong.
Chris
Perry de Valpine
May 29, 2019, 6:08:55 PM5/29/19
to Lokesh Arya, nimble-users
Hi Lokesh,
I'm following up on your case of using nimbleRcall to call dcomp. I'm glad the workaround of coding dcomp as a nimbleFunction worked. I tracked down why the nimbleRcall didn't work and have fixed it on our devel branch. If you need it, you can install by:
library(devtools)
install_github("nimble-dev/nimble", ref = "devel", subdir = "packages/nimble")
install_github("nimble-dev/nimble", ref = "devel", subdir = "packages/nimble")
It will be included in version 0.8.0, which we are working to release soon.
The issue with the nimbleRcall involved management of R's random number generator (RNG). R provides a system for C code to obtain the RNG seed, use RNG calls as needed, and then put the RNG seed back into R when done. Both nimble and Rcpp (via Rcpp::export) automatically do this for all calls into C code. Unfortunately, the C++ generated for a nimbleRcall neglected to manage the RNG seed properly before and after handing control to R. Since in your case the R code called by your nimbleRcall went on to call into C compiled via Rcpp (i.e. dcomp), there was incorrect management of the RNG seed.
-Perry
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/01f86385-1446-460e-8443-18a9ddc4ed3b%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages