[精彩WinBUGS答问]WinBUGS error with zero values in binomial distribution

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

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

############################################

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：distribution Binomial winbugs WINBUG nomial following problems message appears simple

For your 1st case, the problem is when N[i]=0 !

This problem is easy to reproduce and it's straightforward to test his conclusion. E.g., compile this model

model {
    for (i in 1:sites) {
        N[i] ~ dpois(lambda)
        for (j in 1:sample) {
            y[i, j] ~ dbin(p, N[i])
        }
    }
}
with a simple set of data such as

list(
  y = structure(
    .Data = c(0,0,0,0,0,0,0,0),
    .Dim = c(4,2)
  ),
  lambda = 1,
  sites = 4, sample = 2, p = 0.5
)
(Note that all data values are zero.) When you try to generate initial values for the chains, your error is raised:

value of order of binomial y[1,1] must be greater than zero

Now change lambda from a small value to a large one (for which N[i]==0 will essentially never happen):

...
lambda = 100,
...
The model runs fine.

You can also vary the data c(0,0,0,0,0,0,0,0) to confirm that their values do not affect the error. Ergo, the error must be due to either to an invalid value of N[i] or p. Clearly p (= 0.5) is valid, QED.