Skip to content

Commit

Permalink
dNmixture BBNB (non-AD) missing value support
Browse files Browse the repository at this point in the history
  • Loading branch information
@kenkellner
kenkellner committed Sep 3, 2024
1 parent 068cd1d commit 23dc177
Show file tree
Hide file tree
Showing 3 changed files with 168 additions and 10 deletions.
38 changes: 32 additions & 6 deletions R/dNmixture.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -714,9 +714,22 @@ dNmixture_BBNB_v <- nimbleFunction(
if (Nmin == -1 | Nmax == -1) {
stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
}
Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
logProb <- dNmixture_BBNB_steps(x, beta-x,lambda,theta,s,Nmin,Nmax,
dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE))

max_x <- 0
any_not_na <- FALSE
for (i in 1:length(x)){
if(!is.na(x[i])){
if(x[i] > max_x) max_x <- x[i]
any_not_na <- TRUE
}
}
Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x

logProb <- 0
if(any_not_na){
logProb <- dNmixture_BBNB_steps(x, beta-x, lambda, theta, s, Nmin, Nmax,
dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE))
}
if (log) return(logProb)
else return(exp(logProb))
returnType(double())
Expand Down Expand Up @@ -752,9 +765,22 @@ dNmixture_BBNB_s <- nimbleFunction(
if (Nmin == -1 | Nmax == -1) {
stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
}
Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
logProb <- dNmixture_BBNB_steps(x, beta-x,lambda,theta,s,Nmin,Nmax,
dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE))

max_x <- 0
any_not_na <- FALSE
for (i in 1:length(x)){
if(!is.na(x[i])){
if(x[i] > max_x) max_x <- x[i]
any_not_na <- TRUE
}
}
Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x

logProb <- 0
if(any_not_na){
logProb <- dNmixture_BBNB_steps(x, beta-x, lambda, theta, s, Nmin, Nmax,
dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE))
}
if (log) return(logProb)
else return(exp(logProb))
returnType(double())
Expand Down
15 changes: 11 additions & 4 deletions R/utils.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -303,15 +303,22 @@ dNmixture_BBNB_steps <- nimbleFunction(
numN <- NmaxC - NminC + 1 - 1 # remember...
prods <- rep(0, numN)
# N.B. alpha+beta == s
for (i in (NminC + 1):NmaxC)
prods[i - NminC] <- prod(i * (i - 1 + beta_m_x) / ((i - x) * (s + i - 1))) *
((1 - pNB) * (i + r - 1) / i)
for (i in (NminC + 1):NmaxC){
prodi <- 1
for (j in 1:length(x)){
xj <- ADbreak(x[j])
if(!is.na(xj)){
prodi <- prodi * (i * (i - 1 + beta_m_x[j]) / ((i - x[j]) * (s + i - 1)))
}
prods[i - NminC] <- prodi * ((1 - pNB) * (i + r - 1) / i)
}
}
ff <- log(prods)
log_fac <- nimNmixPois_logFac(numN, ff, max_index)
logProb <- logProb + log_fac
}
return(logProb)
returnType(double())
},
buildDerivs = list(run = list(ignore = c("i")))
buildDerivs = list(run = list(ignore = c("i","j","xj")))
)
125 changes: 125 additions & 0 deletions tests/testthat/test-Nmixture.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1528,6 +1528,65 @@ test_that("dNmixture_BBNB_v works",
CMlProbX <- cm$getLogProb("x")
expect_equal(CMlProbX, lProbX)

# Missing values
# Uncompiled calculation
x <- c(1, 0, NA, 3, 0)
Nmin <- max(x, na.rm=TRUE)

probX <- dNmixture_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len)

# Manually calculate the correct answer
alpha <- prob * s
beta <- s - prob * s
r <- 1 / theta
pNB <- 1 / (1 + theta * lambda)

correctProbX <- 0
for (N in Nmin:Nmax) {
correctProbX <- correctProbX + dnbinom(N, size = r, prob = pNB) *
prod(dBetaBinom_v(x, N, shape1 = alpha, shape2 = beta), na.rm=TRUE)
}

expect_equal(probX, correctProbX)

# Uncompiled log probability
lProbX <- dNmixture_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE)
lCorrectProbX <- log(correctProbX)
expect_equal(lProbX, lCorrectProbX)

# Compilation and compiled calculations
CdNmixture_BBNB_v <- compileNimble(dNmixture_BBNB_v)
CprobX <- CdNmixture_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len)
expect_equal(CprobX, probX)

ClProbX <- CdNmixture_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE)
expect_equal(ClProbX, lProbX)

# Use in Nimble model
m <- nimbleModel(code = nc,
data = list(x = x),
inits = list(lambda = lambda,
prob = prob,
s = s, theta = theta),
constants = list(Nmin = Nmin, Nmax = Nmax,
len = len))
m$calculate()
MlProbX <- m$getLogProb("x")
expect_equal(MlProbX, lProbX)

# Compiled model
cm <- compileNimble(m)
cm$calculate()
CMlProbX <- cm$getLogProb("x")
expect_equal(CMlProbX, lProbX)

x <- as.numeric(rep(NA,5))
probX <- dNmixture_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len)
expect_equal(probX, 1)

CprobX <- CdNmixture_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len)
expect_equal(CprobX, 1)

# Test imputing value for all NAs
xNA <- c(NA, NA, NA, NA, NA)
mNA <- nimbleModel(nc, data = list(x = xNA),
Expand Down Expand Up @@ -1670,6 +1729,63 @@ test_that("dNmixture_BBNB_s works",
CMlProbX <- cm$getLogProb("x")
expect_equal(CMlProbX, lProbX)

# Missing values
# Uncompiled calculation
x <- c(1, 0, NA, 3, 0)
Nmin <- max(x, na.rm=TRUE)
probX <- dNmixture_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len)

# Manually calculate the correct answer
alpha <- prob * s
beta <- s - prob * s
r <- 1 / theta
pNB <- 1 / (1 + theta * lambda)

correctProbX <- 0
for (N in Nmin:Nmax) {
correctProbX <- correctProbX + dnbinom(N, size = r, prob = pNB) *
prod(dBetaBinom_s(x, N, shape1 = alpha, shape2 = beta), na.rm=TRUE)
}

expect_equal(probX, correctProbX)

# Uncompiled log probability
lProbX <- dNmixture_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE)
lCorrectProbX <- log(correctProbX)
expect_equal(lProbX, lCorrectProbX)

# Compilation and compiled calculations
CprobX <- CdNmixture_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len)
expect_equal(CprobX, probX)

ClProbX <- CdNmixture_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE)
expect_equal(ClProbX, lProbX)

m <- nimbleModel(code = nc,
data = list(x = x),
inits = list(lambda = lambda,
prob = prob,
s = s, theta = theta),
constants = list(Nmin = Nmin, Nmax = Nmax,
len = len))
m$calculate()
MlProbX <- m$getLogProb("x")
expect_equal(MlProbX, lProbX)

# Compiled model
cm <- compileNimble(m)
cm$calculate()
CMlProbX <- cm$getLogProb("x")
expect_equal(CMlProbX, lProbX)

x <- as.numeric(rep(NA,5))
probX <- dNmixture_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len)
expect_equal(probX, 1)

# Compilation and compiled calculations
CprobX <- CdNmixture_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len)
expect_equal(CprobX, 1)

# Test imputing value for all NAs
xNA <- c(NA, NA, NA, NA, NA)
mNA <- nimbleModel(nc, data = list(x = xNA),
Expand Down Expand Up @@ -1813,6 +1929,15 @@ test_that("dNmixture_BBNB_oneObs works",
CMlProbX <- cm$getLogProb("x")
expect_equal(CMlProbX, lProbX)

# Missing value
# Uncompiled calculation
x <- as.numeric(NA)
probX <- dNmixture_BBNB_oneObs(x, lambda, theta, prob, s, Nmin, Nmax)
expect_equal(probX, 1)

CprobX <- CdNmixture_BBNB_oneObs(x, lambda, theta, prob, s, Nmin, Nmax)
expect_equal(CprobX, 1)

# Test imputing value for all NAs
xNA <- NA
mNA <- nimbleModel(nc, data = list(x = xNA),
Expand Down

0 comments on commit 23dc177

Please sign in to comment.