WinBUGS error with zero values in binomial distribution: value of order of binomial <expr> must be greater than zero
up vote
-1
down vote
favorite
It seems that WinBUGS has problems if it has only zero draws from one binomial distribution:
1. case - simple model
for (i in 1:sites) {
N ~ dpois(lambda)
for (j in 1:sample) {
y[i, j] ~ dbin(p, N)
}
}
If all values y[i,] are zero, then the following error message appears:
value of order of binomial y[53,1] must be greater than zero
EDIT: note that N was not zero in generated data upon failure.
2. case - advanced removal model
for (i in 1:sites) {
M ~ dpois(lambda)
for (j in 1:sample) {
y[i, j, 1] ~ dbin(p, M)
after_removal[i, j] <- M - y[i, j, 1]
y[i, j, 2] ~ dbin(p, after_removal[i, j])
}
}
Here it seems that if some y[i, j, 2] is zero then the error message (same as in case 1) appears!
Note that in this case, I found a workaround: I reparametrized the model to different one which I hope is equivalent:
for (i in 1:sites) {
M ~ dpois(lambda)
for (j in 1:sample) {
y[i, j, 1] ~ dbin(p, M)
obs[i, j] <- y[i, j, 1] + y[i, j, 2]
obs[i, j] ~ dbin(p_komb, M)
}
}
p_komb <- p + (1 - p) * p
... anyway an explanation, clean solution and solution for case 1 wanted:
Why this happens? Is it a WinBUGS bug? How can it be overriden?I have WinBUGS 1.4.3 (August 2007) with immortality patch installed. Below is complete reproducible code for R and package R2WinBUGS (with data generation):
1. case
require(vcd)
sites <- 120 # 60
mean_N <- 16
N <- rpois(sites, mean_N)
p <- 0.4
sample <- 3 # 3
y = matrix(nrow = sites, ncol = sample)
for (i in 1:sites) {
y[i,] = rbinom(sample, N, p)
}
y[20,] = 0
sink("tmp_bugs/model.txt")
cat("
model {
# likelihood
for (i in 1:sites) {
N ~ dpois(lambda)
for (j in 1:sample) {
y[i, j] ~ dbin(p, N)
}
}
# derived parameters
Ntot <- sum(N[])
# priors
p <- 1/(1+exp(-logit_p))
tau <- 1/(4 * 4)
logit_p ~ dnorm(0, tau)
lambda ~ dunif(0, 100)
}
")
sink()
win.data = list(y = y, sample = sample, sites = sites)
#win.data = list(y = y)
#inits = function () { list(N = apply(y, 1, max), p = mean(apply(y, 1, mean)/apply(y, 1, max))) }
inits = function () { list(N = apply(y, 1, max), logit_p = rnorm(1, 0, 4)) }
params = c("N", "p", "Ntot", "lambda")
ni <- 500
nt <- 12
nb <- 200
nc <- 3
out <- bugs(win.data, inits, params, "model.txt",
nc, ni, nb, nt, bugs.directory = bugs.dir,
working.directory = paste(getwd(), "/tmp_bugs/", sep = ""),
debug = TRUE
)
print(out, dig = 3)
par(mfrow = c(2, 2))
hist(out$sims.list$p, breaks = 100)
abline(v = out$mean$p, col = "red", lwd = 2)
abline(v = p, col = "green", lwd = 2)
#lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)
lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4)
legend("topleft", c("estimated", "real", "95% cred. int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7)
#hist(mean(y) / out$sims.list$p, breaks = 1000, xlim = c(0, 50))
hist(out$sims.list$Ntot, breaks = 100)
abline(v = out$mean$Ntot, col = "red", lwd = 2)
abline(v = sum(N), col = "green", lwd = 2)
#lines(quantile(out$sims.list$Ntot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)
lines(quantile(out$sims.list$Ntot, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4)
legend("topright", c("estimated total N", "real total N", "95% credible int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7)
2. case
(one of the y[i, j, 2] will likely be zero; if not, it can be assigned directly)
require(vcd)
sites <- 120 # 60
mean_M <- 16
M <- rpois(sites, mean_M)
p <- 0.4 #0.64
sample <- 2 # 3
y = rep(NA, sites * sample * 2)
dim(y) = c(sites, sample, 2)
for (i in 1:sites) {
# obs[i,] = rbinom(sample, M, p)
for (j in 1:sample) {
y[i,j,1] = rbinom(1, M, p)
y[i,j,2] = rbinom(1, M - y[i,j,1], p)
}
}
y_sample_total = apply(y, c(1, 2), sum)
obs = y_sample_total
############################################
扫码加好友，拉您进群
全部版块
我的主页
收藏
