From 18d46deb0ef1fd7d5339dc20876e6d287c4782b2 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Tue, 13 Aug 2024 14:04:19 -0400 Subject: [PATCH 01/11] Handle NAs in dOcc --- R/dOcc.R | 26 +++++++++++++++++----- tests/testthat/test-AD.R | 47 +++++++++++++++++++++++++++++++++++++++ tests/testthat/test-Occ.R | 46 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 113 insertions(+), 6 deletions(-) diff --git a/R/dOcc.R b/R/dOcc.R index b7d6233..672c34d 100644 --- a/R/dOcc.R +++ b/R/dOcc.R @@ -137,12 +137,19 @@ dOcc_s <- nimbleFunction( log = logical(0, default = 0)) { if (len != 0) if (len != length(x)) stop("Argument 'len' must match length of data, or be 0.") returnType(double(0)) - logProb_x_given_occupied <- sum(dbinom(x, prob = probDetect, size = 1, log = TRUE)) - prob_x_given_unoccupied <- sum(x) == 0 + logProb_x_given_occupied <- 0 + prob_x_given_unoccupied <- 1 + for(i in 1:length(x)) { + xi <- ADbreak(x[i]) + if(!is.na(xi)) { # Handle missing values + logProb_x_given_occupied <- logProb_x_given_occupied + dbinom(x[i], prob = probDetect, size = 1, log = TRUE) + if(xi==1) prob_x_given_unoccupied <- 0 + } + } prob_x <- exp(logProb_x_given_occupied) * probOcc + prob_x_given_unoccupied * (1 - probOcc) if (log) return(log(prob_x)) return(prob_x) - }, buildDerivs = TRUE + }, buildDerivs = list(run = list(ignore = c("i", "xi"))) ) #' @export @@ -156,12 +163,19 @@ dOcc_v <- nimbleFunction( if (len != 0) if (len != length(x)) stop("Argument 'len' must match length of data, or be 0.") if (length(x) != length(probDetect)) stop("Length of data does not match length of detection vector.") returnType(double(0)) - logProb_x_given_occupied <- sum(dbinom(x, prob = probDetect, size = 1, log = TRUE)) - prob_x_given_unoccupied <- sum(x) == 0 + logProb_x_given_occupied <- 0 + prob_x_given_unoccupied <- 1 + for(i in 1:length(x)) { + xi <- ADbreak(x[i]) + if(!is.na(xi)) { # Handle missing values + logProb_x_given_occupied <- logProb_x_given_occupied + dbinom(x[i], prob = probDetect[i], size = 1, log = TRUE) + if(xi==1) prob_x_given_unoccupied <- 0 + } + } prob_x <- exp(logProb_x_given_occupied) * probOcc + prob_x_given_unoccupied * (1 - probOcc) if (log) return(log(prob_x)) return(prob_x) - }, buildDerivs = TRUE + }, buildDerivs = list(run = list(ignore = c("i", "xi"))) ) #' @export diff --git a/tests/testthat/test-AD.R b/tests/testthat/test-AD.R index 3089fc6..dfdf47f 100644 --- a/tests/testthat/test-AD.R +++ b/tests/testthat/test-AD.R @@ -51,6 +51,33 @@ test_that("dOcc works with AD", v1_case1, v2_case1, 0:2) # lots of output numbers with no warning messages means it passes. + # Missing values + dat2 <- c(1,NA,0,0) # A vector of observations + + nc <- nimbleCode({ + x[1:4] ~ dOcc_s(probOcc, probDetect, len = 4) + probOcc ~ dunif(0,1) + probDetect ~ dunif(0,1) + }) + + Rmodel_na <- nimbleModel(nc, data = list(x = dat2), + inits = list(probOcc = probOcc, + probDetect = probDetect), + buildDerivs=TRUE) + + Cmodel_na <- compileNimble(Rmodel_na) + + nodesList_case_na <- + setup_update_and_constant_nodes_for_tests(Rmodel_na, c('probOcc', 'probDetect')) + v1_case1 <- list(arg1 = c(0.6, 0.4)) # taping values for probOcc and probDetect + v2_case1 <- list(arg1 = c(0.65, 0.35)) # testing values for probOcc and probDetect + RCrelTol = c(1e-15, 1e-8, 1e-3, 1e-14) + + res_na <- model_calculate_test_case(Rmodel_na, Cmodel_na, + model_calculate_test, nodesList_case_na, + v1_case1, v2_case1, + 0:2) # lots of output numbers with no warning messages means it passes. + ##################### #### dOcc_v case #### @@ -86,6 +113,26 @@ test_that("dOcc works with AD", model_calculate_test, nodesList_case1, v1_case1, v2_case1, 0:2) + + # Missing values + x2 <- c(1,0,NA,1,0) + Rmodel_na <- nimbleModel(nc, data = list(x = x2), + inits = list(probOcc = probOcc, + probDetect = probDetect), + buildDerivs=TRUE) + Rmodel_na$calculate() + + Cmodel_na <- compileNimble(Rmodel_na) + + nodesList_case_na <- setup_update_and_constant_nodes_for_tests(Rmodel_na, + c('probOcc', Rmodel_na$expandNodeNames('probDetect[1:5]'))) + v1_case1 <- list(arg1 = c(probOcc, probDetect)) # taping values for probOcc and probDetect + v2_case1 <- list(arg1 = c(probOcc2, probDetect2)) # testing values for probOcc and probDetect + + res <- model_calculate_test_case(Rmodel_na, Cmodel_na, + model_calculate_test, nodesList_case_na, + v1_case1, v2_case1, + 0:2) }) test_that ("dNmixture works with AD", { diff --git a/tests/testthat/test-Occ.R b/tests/testthat/test-Occ.R index d570241..48d7e80 100644 --- a/tests/testthat/test-Occ.R +++ b/tests/testthat/test-Occ.R @@ -40,6 +40,29 @@ test_that("dOcc_s and rOcc_s work", { ClProbX <- CdOcc_s(x, probOcc, probDetect, log = TRUE) expect_equal(ClProbX, lProbX) + # Test missing value handling + x3 <- c(1, NA, 0, 1, 0) + probX3 <- dOcc_s(x3, probOcc, probDetect) + correctProbX3 <- + probOcc * probDetect^2 * + (1 - probDetect)^2 + expect_equal(probX3, correctProbX3) + CprobX3 <- CdOcc_s(x3, probOcc, probDetect) + expect_equal(CprobX3, correctProbX3) + + x4 <- c(1, NA, NA, NA, NA) # single-visit + probX4 <- dOcc_s(x4, probOcc, probDetect) + correctProbX4 <- probOcc * probDetect + expect_equal(probX4, correctProbX4) + CprobX4 <- CdOcc_s(x4, probOcc, probDetect) + expect_equal(CprobX4, correctProbX4) + + x5 <- as.numeric(c(NA, NA, NA, NA, NA)) # all missing (as.numeric necessary!) + probX5 <- dOcc_s(x5, probOcc, probDetect) + expect_equal(probX5, 1) + CprobX5 <- CdOcc_s(x5, probOcc, probDetect) + expect_equal(CprobX5, 1) + set.seed(1) nSim <- 10 xSim <- matrix(nrow = nSim, ncol = 5) @@ -143,6 +166,29 @@ test_that("dOcc_v works", { ClProbX <- CdOcc_v(x, probOcc, probDetect, log = TRUE) expect_equal(ClProbX, lProbX) + # Test missing value handling + x3 <- c(1, NA, 0, 1, 0) + probX3 <- dOcc_v(x3, probOcc, probDetect) + correctProbX3 <- + probOcc * prod(probDetect[which(!is.na(x3) & x3 == 1)]) * + prod(1 - probDetect[which(!is.na(x3) & x3 == 0)]) + expect_equal(probX3, correctProbX3) + CprobX3 <- CdOcc_v(x3, probOcc, probDetect) + expect_equal(CprobX3, correctProbX3) + + x4 <- c(1, NA, NA, NA, NA) # single-visit + probX4 <- dOcc_v(x4, probOcc, probDetect) + correctProbX4 <- probOcc[1] * probDetect[1] + expect_equal(probX4, correctProbX4) + CprobX4 <- CdOcc_v(x4, probOcc, probDetect) + expect_equal(CprobX4, correctProbX4) + + x5 <- as.numeric(c(NA, NA, NA, NA, NA)) # all missing (as.numeric necessary!) + probX5 <- dOcc_v(x5, probOcc, probDetect) + expect_equal(probX5, 1) + CprobX5 <- CdOcc_v(x5, probOcc, probDetect) + expect_equal(CprobX5, 1) + set.seed(1) nSim <- 10 xSim <- matrix(nrow = nSim, ncol = 5) From 372043dfdf4b9fbdca9c1606523e1dc37b00c06c Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 23 Aug 2024 08:50:20 -0400 Subject: [PATCH 02/11] Handle NAs in dBetaBinom --- R/dBetaBinom.R | 26 +++++++++++++--------- tests/testthat/test-BetaBinom.R | 38 +++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 10 deletions(-) diff --git a/R/dBetaBinom.R b/R/dBetaBinom.R index 1a24682..0c76153 100644 --- a/R/dBetaBinom.R +++ b/R/dBetaBinom.R @@ -80,17 +80,20 @@ dBetaBinom_v <- nimbleFunction( logprob <- 0 lgNp1 <- lgamma(N+1) for (i in 1:length(x)) { - logprob <- logprob + - nimBetaFun(a = x[i] + shape1[i], b = N - x[i] + shape2[i], log = TRUE) - - nimBetaFun(a = shape1[i], b = shape2[i], log = TRUE) + - lgNp1 - (lgamma(x[i] + 1) + lgamma(N - x[i] + 1)) + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + logprob <- logprob + + nimBetaFun(a = x[i] + shape1[i], b = N - x[i] + shape2[i], log = TRUE) - + nimBetaFun(a = shape1[i], b = shape2[i], log = TRUE) + + lgNp1 - (lgamma(x[i] + 1) + lgamma(N - x[i] + 1)) + } } if (log) return(logprob) return(exp(logprob)) returnType(double(0)) }, - buildDerivs = list(run = list(ignore = 'i')) + buildDerivs = list(run = list(ignore = c('i','xi'))) ) #' @rdname dBetaBinom @@ -106,16 +109,19 @@ dBetaBinom_s <- nimbleFunction( lgNp1 <- lgamma(N+1) lbs1s2 <- nimBetaFun(a = shape1, b = shape2, log = TRUE) for (i in 1:length(x)) { - logprob <- logprob + - nimBetaFun(a = x[i] + shape1, b = N - x[i] + shape2, log = TRUE) - - lbs1s2 + - lgNp1 - (lgamma(x[i] + 1) + lgamma(N - x[i] + 1)) + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + logprob <- logprob + + nimBetaFun(a = x[i] + shape1, b = N - x[i] + shape2, log = TRUE) - + lbs1s2 + + lgNp1 - (lgamma(x[i] + 1) + lgamma(N - x[i] + 1)) + } } if (log) return(logprob) return(exp(logprob)) returnType(double(0)) }, - buildDerivs = list(run=list(ignore = 'i')) + buildDerivs = list(run=list(ignore = c('i','xi'))) ) #' @rdname dBetaBinom diff --git a/tests/testthat/test-BetaBinom.R b/tests/testthat/test-BetaBinom.R index 8ecf811..b80b44e 100644 --- a/tests/testthat/test-BetaBinom.R +++ b/tests/testthat/test-BetaBinom.R @@ -67,6 +67,25 @@ test_that("dBetaBinom_v works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + xna <- c(4, NA, 8, 0, 3) + probXna <- dBetaBinom_v(xna, N, shape1, shape2) + len <- 5 + + correctProbXna <- prod( + choose(N, xna) * beta(xna + shape1, N - xna + shape2) / + beta(shape1, shape2) + , na.rm=TRUE) + expect_equal(probXna, correctProbXna) + + CprobXna <- CdBetaBinom_v(xna, N, shape1, shape2, len = len) + expect_equal(CprobXna, probXna) + + # All NAs + xna <- as.numeric(rep(NA, 5)) + expect_equal(CdBetaBinom_v(xna, N, shape1, shape2, len=len), 1) + expect_equal(CdBetaBinom_v(xna, N, shape1, shape2, len=len, log=TRUE), 0) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -179,6 +198,25 @@ test_that("dBetaBinom_s works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + xna <- c(4, NA, 8, 0, 3) + probXna <- dBetaBinom_s(xna, N, shape1, shape2) + len <- 5 + + correctProbXna <- prod( + choose(N, xna) * beta(xna + shape1, N - xna + shape2) / + beta(shape1, shape2) + , na.rm=TRUE) + expect_equal(probXna, correctProbXna) + + CprobXna <- CdBetaBinom_s(xna, N, shape1, shape2, len = len) + expect_equal(CprobXna, probXna) + + # All NAs + xna <- as.numeric(rep(NA, 5)) + expect_equal(CdBetaBinom_s(xna, N, shape1, shape2, len=len), 1) + expect_equal(CdBetaBinom_s(xna, N, shape1, shape2, len=len, log=TRUE), 0) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), From efbcfd74625b8cbca9caf4ae5c7c4d2b7106eab4 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 23 Aug 2024 15:36:32 -0400 Subject: [PATCH 03/11] dNmixture / dNmixture-AD na handling --- R/dNmixture.R | 56 +++++++++++++-- R/dNmixtureAD.R | 67 +++++++++++++++--- R/utils.R | 11 ++- tests/testthat/test-AD.R | 36 ++++++++++ tests/testthat/test-Nmixture.R | 86 ++++++++++++++++++++++++ tests/testthat/test-NmixtureADnoDerivs.R | 50 ++++++++++++++ 6 files changed, 288 insertions(+), 18 deletions(-) diff --git a/R/dNmixture.R b/R/dNmixture.R index 4edaa23..72e3022 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -245,9 +245,31 @@ dNmixture_v <- nimbleFunction( Nmin <- min(x + qpois(0.00001, lambda * (1 - prob))) if (Nmax == -1) Nmax <- max(x + qpois(0.99999, lambda * (1 - prob))) - Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x - logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, sum(log(1-prob)), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE))) + + max_x <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + if(x[i] > max_x) max_x <- x[i] + } + } + Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + if(!is.na(x[i])){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob[i]) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob[i], log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, sum_log_one_m_prob, + sum_log_dbinom) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) @@ -277,9 +299,31 @@ dNmixture_s <- nimbleFunction( Nmin <- min(x + qpois(0.00001, lambda * (1 - prob))) if (Nmax == -1) Nmax <- max(x + qpois(0.99999, lambda * (1 - prob))) - Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x - logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, log(1-prob)*len, - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE))) + + max_x <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + if(x[i] > max_x) max_x <- x[i] + } + } + Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + if(!is.na(x[i])){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob, log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, sum_log_one_m_prob, + sum_log_dbinom) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) diff --git a/R/dNmixtureAD.R b/R/dNmixtureAD.R index 76e13d8..2347842 100644 --- a/R/dNmixtureAD.R +++ b/R/dNmixtureAD.R @@ -77,14 +77,38 @@ dNmixtureAD_v <- nimbleFunction( if (log) return(-Inf) else return(0) if ((Nmin == -1) | (Nmax == -1)) stop("Must provide Nmin and Nmax in AD version of dNmixture distributions") - Nmin <- ADbreak(max( max(x), Nmin )) - logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, sum(log(1-prob)), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE)), - usingAD=TRUE) + + max_x <- 0 + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + } + } + Nmin <- ADbreak(max( max_x, Nmin )) + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + # You can't combine this with the loop above because we need to know the final Nmin + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob[i]) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob[i], log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, sum_log_one_m_prob, + sum_log_dbinom, usingAD=TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs = list(run = list()) + }, buildDerivs = list(run = list(ignore=c("i","xi"))) ) #' @rdname dNmixtureAD @@ -102,14 +126,37 @@ dNmixtureAD_s <- nimbleFunction( if (log) return(-Inf) else return(0) if ((Nmin == -1) | (Nmax == -1)) stop("Must provide Nmin and Nmax in AD version of dNmixture distributions") - Nmin <- ADbreak(max( max(x), Nmin )) - logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, len*log(1-prob), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE)), - usingAD=TRUE) + + max_x <- 0 + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + } + } + Nmin <- ADbreak(max( max_x, Nmin )) + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob, log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_steps(x, lambda, Nmin, Nmax, sum_log_one_m_prob, + sum_log_dbinom, usingAD=TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs = list(run = list()) + }, buildDerivs = list(run = list(ignore=c("i","xi"))) ) NULL diff --git a/R/utils.R b/R/utils.R index ef511ee..4aab8a5 100644 --- a/R/utils.R +++ b/R/utils.R @@ -135,7 +135,14 @@ dNmixture_steps <- nimbleFunction( numN <- NmaxC - NminC + 1 - 1 ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1 prods <- rep(0, numN) for (i in (NminC + 1):NmaxC) { - prods[i - NminC] <- prod(i/(i - x)) / i + prodi <- 1 + for (j in 1:length(x)){ + xj <- ADbreak(x[j]) + if(!is.na(xj)){ + prodi <- prodi * (i / (i - x[j])) + } + } + prods[i - NminC] <- prodi / i } ff <- log(lambda) + sum_log_one_m_prob + log(prods) log_fac <- nimNmixPois_logFac(numN, ff, max_index) @@ -144,7 +151,7 @@ dNmixture_steps <- nimbleFunction( return(logProb) returnType(double()) }, - buildDerivs = list(run = list(ignore = c("i"))) + buildDerivs = list(run = list(ignore = c("i","j","xj"))) ) ##### N-mixture extensions ##### diff --git a/tests/testthat/test-AD.R b/tests/testthat/test-AD.R index dfdf47f..781eb20 100644 --- a/tests/testthat/test-AD.R +++ b/tests/testthat/test-AD.R @@ -171,6 +171,24 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # Missing values + xna <- c(7, 7, NA, 9, 10) + Rmodel_na <- nimbleModel(nc, data = list(x = xna), + inits = list(prob = prob, + lambda = lambda), + buildDerivs=TRUE) + Rmodel_na$calculate() + + Cmodel_na <- compileNimble(Rmodel_na) + Cmodel_na$calculate() + + nodesList_case1_na <- setup_update_and_constant_nodes_for_tests(Rmodel_na, c('prob', 'lambda')) + + res_na <- model_calculate_test_case(Rmodel_na, Cmodel_na, + model_calculate_test, nodesList_case1_na, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ########################## #### dNmixture_BNB_s case #### @@ -325,6 +343,24 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # Missing values + xna <- c(7, 7, NA, 9, 10) + Rmodel_na <- nimbleModel(nc, data = list(x = xna), + inits = list(prob = prob, + lambda = lambda), + buildDerivs=TRUE) + Rmodel_na$calculate() + + Cmodel_na <- compileNimble(Rmodel_na) + Cmodel_na$calculate() + + nodesList_case1_na <- setup_update_and_constant_nodes_for_tests(Rmodel_na, c('prob', 'lambda')) + + res_na <- model_calculate_test_case(Rmodel_na, Cmodel_na, + model_calculate_test, nodesList_case1_na, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ############################## #### dNmixture_BNB_v case #### diff --git a/tests/testthat/test-Nmixture.R b/tests/testthat/test-Nmixture.R index 32f1e81..c4da0a7 100644 --- a/tests/testthat/test-Nmixture.R +++ b/tests/testthat/test-Nmixture.R @@ -94,6 +94,49 @@ test_that("dNmixture_v works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing value handling + xna <- c(1, 0, NA, 3, 0) + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + probXna <- dNmixture_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, correctProbXna) + CprobXna <- CdNmixture_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, CprobXna) + + m_na <- nimbleModel(code = nc, + data = list(x = xna), + inits = list(lambda = lambda, + prob = prob), + constants = list(Nmin = Nmin, Nmax = Nmax, + len = len)) + m_na$calculate() + MlProbX_na <- m_na$getLogProb("x") + expect_equal(MlProbX_na, log(correctProbXna)) + + cm_na <- compileNimble(m_na) + cm_na$calculate() + CMlProbX_na <- cm_na$getLogProb("x") + expect_equal(CMlProbX_na, log(correctProbXna)) + + xna <- c(1, NA, NA, NA, NA) + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + probXna <- dNmixture_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, correctProbXna) + CprobXna <- CdNmixture_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, CprobXna) + + xna <- as.numeric(c(NA, NA, NA, NA, NA)) + probXna <- dNmixture_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, 1) + CprobXna <- CdNmixture_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, CprobXna) + expect_equal(CdNmixture_v(xna, lambda, prob, Nmin, Nmax, len, log=TRUE), 0) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -217,6 +260,49 @@ test_that("dNmixture_s works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing value handling + xna <- c(1, 0, NA, 3, 2) + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + probXna <- dNmixture_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, correctProbXna) + CprobXna <- CdNmixture_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, CprobXna) + + m_na <- nimbleModel(code = nc, + data = list(x = xna), + inits = list(lambda = lambda, + prob = prob), + constants = list(Nmin = Nmin, Nmax = Nmax, + len = len)) + m_na$calculate() + MlProbX_na <- m_na$getLogProb("x") + expect_equal(MlProbX_na, log(correctProbXna)) + + cm_na <- compileNimble(m_na) + cm_na$calculate() + CMlProbX_na <- cm_na$getLogProb("x") + expect_equal(CMlProbX_na, log(correctProbXna)) + + xna <- c(1, NA, NA, NA, NA) + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + probXna <- dNmixture_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, correctProbXna) + CprobXna <- CdNmixture_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, CprobXna) + + xna <- as.numeric(c(NA, NA, NA, NA, NA)) + probXna <- dNmixture_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, 1) + CprobXna <- CdNmixture_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, CprobXna) + expect_equal(CdNmixture_s(xna, lambda, prob, Nmin, Nmax, len, log=TRUE), 0) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), diff --git a/tests/testthat/test-NmixtureADnoDerivs.R b/tests/testthat/test-NmixtureADnoDerivs.R index 67fe741..50d4231 100644 --- a/tests/testthat/test-NmixtureADnoDerivs.R +++ b/tests/testthat/test-NmixtureADnoDerivs.R @@ -87,6 +87,31 @@ test_that("dNmixtureAD_v works uncompiled", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + xna <- c(1, 0, NA, 3, 0) + probXna <- dNmixtureAD_v(xna, lambda, prob, Nmin, Nmax, len) + # Manually calculate the correct answer + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + expect_equal(probXna, correctProbXna) + + CprobXna <- CdNmixtureAD_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(CprobXna, probXna) + + xna <- c(1,NA,NA,NA,NA) + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + probXna <- CdNmixtureAD_v(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, correctProbXna) + + xna <- as.numeric(rep(NA, 5)) + expect_equal(CdNmixtureAD_v(xna, lambda, prob, Nmin, Nmax, len), 1) + expect_equal(CdNmixtureAD_v(xna, lambda, prob, Nmin, Nmax, len, log=TRUE), 0) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -231,6 +256,31 @@ test_that("dNmixtureAD_s works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + xna <- c(1, 0, NA, 3, 2) + probXna <- dNmixtureAD_s(xna, lambda, prob, Nmin, Nmax, len) + # Manually calculate the correct answer + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + expect_equal(probXna, correctProbXna) + + CprobXna <- CdNmixtureAD_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(CprobXna, probXna) + + xna <- c(1,NA,NA,NA,NA) + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dpois(N, lambda) * prod(dbinom(xna, N, prob), na.rm=TRUE) + } + probXna <- CdNmixtureAD_s(xna, lambda, prob, Nmin, Nmax, len) + expect_equal(probXna, correctProbXna) + + xna <- as.numeric(rep(NA, 5)) + expect_equal(CdNmixtureAD_s(xna, lambda, prob, Nmin, Nmax, len), 1) + expect_equal(CdNmixtureAD_s(xna, lambda, prob, Nmin, Nmax, len, log=TRUE), 0) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), From fe20fcd09ce70635e946b48eefa12b296ca95189 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Thu, 29 Aug 2024 14:12:36 -0400 Subject: [PATCH 04/11] dNmixture_BNB / dNmixture_BNB_AD handling --- R/dNmixture.R | 106 +++++++-- R/dNmixtureAD.R | 66 +++++- R/utils.R | 14 +- tests/testthat/test-AD.R | 49 +++++ tests/testthat/test-Nmixture.R | 113 ++++++++++ tests/testthat/test-NmixtureADnoDerivs.R | 267 +++++++++++++++++++++++ 6 files changed, 591 insertions(+), 24 deletions(-) diff --git a/R/dNmixture.R b/R/dNmixture.R index 72e3022..ee2aea6 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -401,13 +401,52 @@ dNmixture_BNB_v <- nimbleFunction( lambda_cond <- omega / (theta_cond * (1 - omega)) r_cond <- 1 / theta_cond pNB_cond <- 1 / (1 + theta_cond * lambda_cond) - if (Nmin == -1) - Nmin <- min(x + qnbinom(0.00001, size = r_cond, prob = pNB_cond)) - if (Nmax == -1) - Nmax <- max(x + qnbinom(0.99999, size = r_cond, prob = pNB_cond)) - Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x - logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax,sum(log(1-prob)), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE))) + + if (Nmax == -1){ + Nmax <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmax_cand <- x[i] + qnbinom(0.99999, size = r_cond[i], prob = pNB_cond[i]) + if(Nmax_cand > Nmax) Nmax <- Nmax_cand + } + } + } + + if(Nmin == -1){ + Nmin <- Nmax + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmin_cand <- x[i] + qnbinom(0.00001, size = r_cond[i], prob = pNB_cond[i]) + if(Nmin_cand < Nmin) Nmin <- Nmin_cand + } + } + } + + max_x <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + if(x[i] > max_x) max_x <- x[i] + } + } + + Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + if(!is.na(x[i])){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob[i]) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob[i], log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax, sum_log_one_m_prob, + sum_log_dbinom) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) @@ -439,13 +478,52 @@ dNmixture_BNB_s <- nimbleFunction( lambda_cond <- omega / (theta_cond * (1 - omega)) r_cond <- 1 / theta_cond pNB_cond <- 1 / (1 + theta_cond * lambda_cond) - if (Nmin == -1) - Nmin <- min(x + qnbinom(0.00001, size = r_cond, prob = pNB_cond)) - if (Nmax == -1) - Nmax <- max(x + qnbinom(0.99999, size = r_cond, prob = pNB_cond)) - Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x - logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax,len*log(1-prob), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE))) + + if (Nmax == -1){ + Nmax <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmax_cand <- x[i] + qnbinom(0.99999, size = r_cond[i], prob = pNB_cond[i]) + if(Nmax_cand > Nmax) Nmax <- Nmax_cand + } + } + } + + if(Nmin == -1){ + Nmin <- Nmax + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmin_cand <- x[i] + qnbinom(0.00001, size = r_cond[i], prob = pNB_cond[i]) + if(Nmin_cand < Nmin) Nmin <- Nmin_cand + } + } + } + + max_x <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + if(x[i] > max_x) max_x <- x[i] + } + } + + Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + if(!is.na(x[i])){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob, log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax, sum_log_one_m_prob, + sum_log_dbinom) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) diff --git a/R/dNmixtureAD.R b/R/dNmixtureAD.R index 2347842..0ae5f09 100644 --- a/R/dNmixtureAD.R +++ b/R/dNmixtureAD.R @@ -210,13 +210,38 @@ dNmixtureAD_BNB_v <- nimbleFunction( if (log) return(-Inf) else return(0) if ((Nmin == -1) | (Nmax == -1)) stop("Must provide Nmin and Nmax in AD version of dNmixture distributions") - Nmin <- ADbreak(max( max(x), Nmin )) - logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax,sum(log(1-prob)), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE))) + + max_x <- 0 + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + } + } + + Nmin <- ADbreak(max( max_x, Nmin )) + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob[i]) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob[i], log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax, sum_log_one_m_prob, + sum_log_dbinom, usingAD = TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs = list(run = list()) + }, buildDerivs = list(run = list(ignore=c("i", "xi"))) ) #' @rdname dNmixtureAD @@ -237,13 +262,38 @@ dNmixtureAD_BNB_s <- nimbleFunction( if (log) return(-Inf) else return(0) if ((Nmin == -1) | (Nmax == -1)) stop("Must provide Nmin and Nmax in AD version of dNmixture distributions") - Nmin <- ADbreak(max( max(x), Nmin )) - logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax,len*log(1-prob), - sum(dbinom(x, size = Nmin, prob = prob, log = TRUE))) + + max_x <- 0 + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + } + } + + Nmin <- ADbreak(max( max_x, Nmin )) + + sum_log_dbinom <- 0 + sum_log_one_m_prob <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + sum_log_one_m_prob <- sum_log_one_m_prob + log(1 - prob) + sum_log_dbinom <- sum_log_dbinom + dbinom(x[i], size = Nmin, prob = prob, log=TRUE) + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BNB_steps(x,lambda,theta,Nmin,Nmax, sum_log_one_m_prob, + sum_log_dbinom, usingAD = TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs = list(run = list()) + }, buildDerivs = list(run = list(ignore=c("i", "xi"))) ) #' @rdname dNmixtureAD diff --git a/R/utils.R b/R/utils.R index 4aab8a5..5747a2e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -190,8 +190,18 @@ dNmixture_BNB_steps <- nimbleFunction( numN <- 0L numN <- NmaxC - NminC + 1 - 1 # remember... prods <- rep(0, numN) - for (i in (NminC + 1):NmaxC) + for (i in (NminC + 1):NmaxC){ prods[i - NminC] <- (i + r - 1) * prod(i/(i - x)) / i + + prodi <- i + r - 1 + for (j in 1:length(x)){ + xj <- ADbreak(x[j]) + if(!is.na(xj)){ + prodi <- prodi * (i / (i - x[j])) + } + } + prods[i - NminC] <- prodi / i + } ff <- log(1 - pNB) + sum_log_one_m_prob + log(prods) log_fac <- nimNmixPois_logFac(numN, ff, max_index) logProb <- logProb + log_fac @@ -199,7 +209,7 @@ dNmixture_BNB_steps <- nimbleFunction( return(logProb) returnType(double()) }, - buildDerivs = list(run = list(ignore = c("i"))) + buildDerivs = list(run = list(ignore = c("i", "j", "xj"))) ) diff --git a/tests/testthat/test-AD.R b/tests/testthat/test-AD.R index 781eb20..6b7cb18 100644 --- a/tests/testthat/test-AD.R +++ b/tests/testthat/test-AD.R @@ -227,6 +227,34 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # NA handling + x <- c(7, 7, NA, 9, 10) + nc <- nimbleCode({ + x[1:5] ~ dNmixtureAD_BNB_s(lambda, prob, theta = theta, + Nmin = 0, Nmax = 100, len = 5) + prob ~ dunif(0, 1) + lambda ~ dunif(0, 100) + }) + + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, + lambda = lambda, + theta = theta), + buildDerivs=TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta')) + v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ############################## #### dNmixture_BBP_s case #### @@ -403,6 +431,27 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # NA handling + x <- c(7, 7, NA, 9, 10) + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, + lambda = lambda, + theta = theta), + buildDerivs=TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta')) + v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ############################## #### dNmixture_BBP_v case #### diff --git a/tests/testthat/test-Nmixture.R b/tests/testthat/test-Nmixture.R index c4da0a7..4c3a0e6 100644 --- a/tests/testthat/test-Nmixture.R +++ b/tests/testthat/test-Nmixture.R @@ -458,6 +458,57 @@ test_that("dNmixture_BNB_v works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + xna <- c(1, 0, NA, 3, 0) + probXna <- dNmixture_BNB_v(xna, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(xna, N, prob), na.rm=TRUE) + } + expect_equal(probXna, correctProbXna) + + # Check compiled version + CprobXna <- CdNmixture_BNB_v(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobXna, correctProbXna) + + m_na <- nimbleModel(code = nc, + data = list(x = xna), + inits = list(lambda = lambda, + prob = prob, + theta = theta), + constants = list(Nmin = Nmin, Nmax = Nmax, + len = len)) + m_na$calculate() + MlProbXna <- m_na$getLogProb("x") + expect_equal(MlProbXna, log(correctProbXna)) + + # Compiled model + cm_na <- compileNimble(m_na) + cm_na$calculate() + CMlProbXna <- cm_na$getLogProb("x") + expect_equal(CMlProbXna, log(correctProbXna)) + + xna <- c(1, NA, NA, NA, NA) + probXna <- dNmixture_BNB_v(xna, lambda, theta, prob, Nmin, Nmax, len) + # Manually calculate the correct answer + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(xna, N, prob), na.rm=TRUE) + } + expect_equal(probXna, correctProbXna) + CprobXna <- CdNmixture_BNB_v(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobXna, correctProbXna) + + xna <- as.numeric(rep(NA, 5)) + probXna <- dNmixture_BNB_v(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(probXna, 1) + CprobXna <- CdNmixture_BNB_v(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobXna, 1) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -592,6 +643,57 @@ test_that("dNmixture_BNB_s works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Test NA handling + xna <- c(1, 0, NA, 3, 0) + probXna <- dNmixture_BNB_s(xna, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(xna, N, prob), na.rm=TRUE) + } + expect_equal(probXna, correctProbXna) + + # Check compiled version + CprobXna <- CdNmixture_BNB_s(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobXna, correctProbXna) + + m_na <- nimbleModel(code = nc, + data = list(x = xna), + inits = list(lambda = lambda, + prob = prob, + theta = theta), + constants = list(Nmin = Nmin, Nmax = Nmax, + len = len)) + m_na$calculate() + MlProbXna <- m_na$getLogProb("x") + expect_equal(MlProbXna, log(correctProbXna)) + + # Compiled model + cm_na <- compileNimble(m_na) + cm_na$calculate() + CMlProbXna <- cm_na$getLogProb("x") + expect_equal(CMlProbXna, log(correctProbXna)) + + xna <- c(1, NA, NA, NA, NA) + probXna <- dNmixture_BNB_s(xna, lambda, theta, prob, Nmin, Nmax, len) + # Manually calculate the correct answer + correctProbXna <- 0 + for (N in Nmin:Nmax) { + correctProbXna <- correctProbXna + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(xna, N, prob), na.rm=TRUE) + } + expect_equal(probXna, correctProbXna) + CprobXna <- CdNmixture_BNB_s(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobXna, correctProbXna) + + xna <- as.numeric(rep(NA, 5)) + probXna <- dNmixture_BNB_s(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(probXna, 1) + CprobXna <- CdNmixture_BNB_s(xna, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobXna, 1) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -725,6 +827,17 @@ test_that("dNmixture_BNB_oneObs works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Test NA handling + xna <- as.numeric(c(NA)) + probXna <- dNmixture_BNB_oneObs(xna, lambda, theta, prob, Nmin, Nmax) + + # Manually calculate the correct answer + expect_equal(probXna, 1) + + # Check compiled version + CprobXna <- CdNmixture_BNB_oneObs(xna, lambda, theta, prob, Nmin, Nmax) + expect_equal(CprobXna, 1) + # Test imputing value for all NAs xNA <- NA mNA <- nimbleModel(nc, data = list(x = xNA), diff --git a/tests/testthat/test-NmixtureADnoDerivs.R b/tests/testthat/test-NmixtureADnoDerivs.R index 50d4231..c32a3da 100644 --- a/tests/testthat/test-NmixtureADnoDerivs.R +++ b/tests/testthat/test-NmixtureADnoDerivs.R @@ -438,6 +438,131 @@ test_that("dNmixtureAD_BNB_v works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + x <- c(1, 0, NA, 3, 0) + probX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(x, N, prob), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Dynamic Nmin/Nmax + expect_error(dynProbX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin = -1, Nmax = -1, len)) + + # Other Nmin / Nmax + Nmin <- 3 + for(Nmax in 3:6) { + dynProbX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin = Nmin, Nmax = Nmax, len) + dynCorrectProbX <- 0 + for (N in Nmin:Nmax) { + dynCorrectProbX <- dynCorrectProbX + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(x, N, prob), na.rm=TRUE) + } + expect_equal(dynProbX, dynCorrectProbX) + } + Nmin <- 0 + Nmax <- 250 + # Compilation and compiled calculations + call_dNmixtureAD_BNB_v <- nimbleFunction( + name = "t3", + run=function(x=double(1), + lambda=double(), + theta=double(), + prob=double(1), + Nmin = double(0), + Nmax = double(0), + len=double(), + log = integer(0, default=0)) { + return(dNmixtureAD_BNB_v(x,lambda,theta,prob,Nmin,Nmax,len,log)) + returnType(double()) + } + ) + + # Compilation and compiled calculations + CdNmixtureAD_BNB_v <- compileNimble(call_dNmixtureAD_BNB_v) + CprobX <- CdNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + probX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lprobX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + expect_equal(ClProbX, lProbX) + + expect_error(CdynProbX <- CdNmixture_BNB_v(x, lambda, theta, prob, Nmin = -1, + Nmax = -1, len)) + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixtureAD_BNB_v(lambda = lambda, prob = prob[1:5], + theta = theta, + Nmin = Nmin, Nmax = Nmax, len = len) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(lambda = lambda, + prob = prob, + 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 <- c(1, NA, NA, NA, NA) + probX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(x, N, prob), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + CprobX <- CdNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + probX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + x <- as.numeric(c(NA, NA, NA, NA, NA)) + probX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbX <- 1 + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + CprobX <- CdNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + probX <- dNmixtureAD_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobX, probX) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -593,6 +718,131 @@ test_that("dNmixtureAD_BNB_s works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + x <- c(1, 0, NA, 3, 0) + probX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(x, N, prob), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Dynamic Nmin/Nmax + expect_error(dynProbX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin = -1, Nmax = -1, len)) + + # Other Nmin / Nmax + Nmin <- 3 + for(Nmax in 3:6) { + dynProbX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin = Nmin, Nmax = Nmax, len) + dynCorrectProbX <- 0 + for (N in Nmin:Nmax) { + dynCorrectProbX <- dynCorrectProbX + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(x, N, prob), na.rm=TRUE) + } + expect_equal(dynProbX, dynCorrectProbX) + } + Nmin <- 0 + Nmax <- 250 + # Compilation and compiled calculations + call_dNmixtureAD_BNB_s <- nimbleFunction( + name = "t3", + run=function(x=double(1), + lambda=double(), + theta=double(), + prob=double(), + Nmin = double(0), + Nmax = double(0), + len=double(), + log = integer(0, default=0)) { + return(dNmixtureAD_BNB_s(x,lambda,theta,prob,Nmin,Nmax,len,log)) + returnType(double()) + } + ) + + # Compilation and compiled calculations + CdNmixtureAD_BNB_s <- compileNimble(call_dNmixtureAD_BNB_s) + CprobX <- CdNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + probX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lprobX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + expect_equal(ClProbX, lProbX) + + expect_error(CdynProbX <- CdNmixture_BNB_s(x, lambda, theta, prob, Nmin = -1, + Nmax = -1, len)) + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixtureAD_BNB_s(lambda = lambda, prob = prob, + theta = theta, + Nmin = Nmin, Nmax = Nmax, len = len) + + }) + + m <- nimbleModel(code = nc, + data = list(x = x), + inits = list(lambda = lambda, + prob = prob, + 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 <- c(1, NA, NA, NA, NA) + probX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dnbinom(N, size = r, prob = pNB) * + prod(dbinom(x, N, prob), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + CprobX <- CdNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + probX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + x <- as.numeric(c(NA, NA, NA, NA, NA)) + probX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + + # Manually calculate the correct answer + correctProbX <- 1 + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + CprobX <- CdNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + probX <- dNmixtureAD_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len) + expect_equal(CprobX, probX) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -742,6 +992,23 @@ test_that("dNmixtureAD_BNB_oneObs works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + x <- as.numeric(NA) + probX <- dNmixtureAD_BNB_oneObs(x, lambda, theta, prob, Nmin, Nmax) + expect_equal(probX, 1) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BNB_oneObs(x, lambda, theta, prob, Nmin, Nmax, log = TRUE) + lCorrectProbX <- log(1) + expect_equal(lProbX, lCorrectProbX) + + # Compilation and compiled calculations + CprobX <- CdNmixtureAD_BNB_oneObs(x, lambda, theta, prob, Nmin, Nmax) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_BNB_oneObs(x, lambda, theta, prob, Nmin, Nmax, log = TRUE) + expect_equal(ClProbX, lProbX) + # Test imputing value for all NAs xNA <- NA mNA <- nimbleModel(nc, data = list(x = xNA), From 4ef1318c1376f2397154af329251f04b0d2ae6f1 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 30 Aug 2024 10:42:20 -0400 Subject: [PATCH 05/11] Handle Nmin/Nmax = -1 in dNmixture --- R/dNmixture.R | 46 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/R/dNmixture.R b/R/dNmixture.R index ee2aea6..9e14f10 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -241,10 +241,25 @@ dNmixture_v <- nimbleFunction( if (log) return(-Inf) else return(0) ## For each x, the conditional distribution of (N - x | x) is pois(lambda * (1-p)) ## We determine the lowest N and highest N at extreme quantiles and sum over those. - if (Nmin == -1) - Nmin <- min(x + qpois(0.00001, lambda * (1 - prob))) - if (Nmax == -1) - Nmax <- max(x + qpois(0.99999, lambda * (1 - prob))) + if (Nmax == -1){ + Nmax <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmax_cand <- x[i] + qpois(0.99999, lambda * (1-prob[i])) + if(Nmax_cand > Nmax) Nmax <- Nmax_cand + } + } + } + + if(Nmin == -1){ + Nmin <- Nmax + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmin_cand <- x[i] + qpois(0.00001, lambda * (1-prob[i])) + if(Nmin_cand < Nmin) Nmin <- Nmin_cand + } + } + } max_x <- 0 for (i in 1:length(x)){ @@ -295,10 +310,25 @@ dNmixture_s <- nimbleFunction( if (log) return(-Inf) else return(0) ## For each x, the conditional distribution of (N - x | x) is pois(lambda * (1-p)) ## We determine the lowest N and highest N at extreme quantiles and sum over those. - if (Nmin == -1) - Nmin <- min(x + qpois(0.00001, lambda * (1 - prob))) - if (Nmax == -1) - Nmax <- max(x + qpois(0.99999, lambda * (1 - prob))) + if (Nmax == -1){ + Nmax <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmax_cand <- x[i] + qpois(0.99999, lambda * (1-prob)) + if(Nmax_cand > Nmax) Nmax <- Nmax_cand + } + } + } + + if(Nmin == -1){ + Nmin <- Nmax + for (i in 1:length(x)){ + if(!is.na(x[i])){ + Nmin_cand <- x[i] + qpois(0.00001, lambda * (1-prob)) + if(Nmin_cand < Nmin) Nmin <- Nmin_cand + } + } + } max_x <- 0 for (i in 1:length(x)){ From 04a6cb2b5b7ca95bb0c01c74032d7bbb81c1fb0c Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 30 Aug 2024 11:11:08 -0400 Subject: [PATCH 06/11] dNmixture_BBP (non-AD) missing value support --- R/dNmixture.R | 47 +++++++++++-- R/utils.R | 14 +++- tests/testthat/test-Nmixture.R | 120 +++++++++++++++++++++++++++++++++ 3 files changed, 172 insertions(+), 9 deletions(-) diff --git a/R/dNmixture.R b/R/dNmixture.R index 9e14f10..6659aff 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -602,9 +602,26 @@ dNmixture_BBP_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_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, - dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE)) + max_x <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + if(x[i] > max_x) max_x <- x[i] + } + } + Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x + + any_not_na <- FALSE + for (i in 1:length(x)){ + if(!is.na(x[i])){ + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, + dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE)) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) @@ -636,9 +653,27 @@ dNmixture_BBP_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_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, - dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE)) + + max_x <- 0 + for (i in 1:length(x)){ + if(!is.na(x[i])){ + if(x[i] > max_x) max_x <- x[i] + } + } + Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x + + any_not_na <- FALSE + for (i in 1:length(x)){ + if(!is.na(x[i])){ + any_not_na <- TRUE + } + } + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, + dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE)) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) diff --git a/R/utils.R b/R/utils.R index 5747a2e..920f96f 100644 --- a/R/utils.R +++ b/R/utils.R @@ -247,8 +247,16 @@ dNmixture_BBP_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))) * (lambda / 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 * (lambda / i) + } + } ff <- log(prods) log_fac <- nimNmixPois_logFac(numN, ff, max_index) logProb <- logProb + log_fac @@ -256,7 +264,7 @@ dNmixture_BBP_steps <- nimbleFunction( return(logProb) returnType(double()) }, - buildDerivs = list(run = list(ignore = c("i"))) + buildDerivs = list(run = list(ignore = c("i","j","xj"))) ) diff --git a/tests/testthat/test-Nmixture.R b/tests/testthat/test-Nmixture.R index 4c3a0e6..4e467e2 100644 --- a/tests/testthat/test-Nmixture.R +++ b/tests/testthat/test-Nmixture.R @@ -981,6 +981,61 @@ test_that("dNmixture_BBP_v works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + # Uncompiled calculation + x <- c(1, 0, NA, 3, 0) + probX <- dNmixture_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) + + # Manually calculate the correct answer + alpha <- prob * s + beta <- s - prob * s + + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dpois(N, lambda) * + prod(dBetaBinom_v(x, N, shape1 = alpha, shape2 = beta), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixture_BBP_v(x, lambda, prob, s, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Compilation and compiled calculations + CdNmixture_BBP_v <- compileNimble(dNmixture_BBP_v) + CprobX <- CdNmixture_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixture_BBP_v(x, lambda, 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), + 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_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixture_BBP_s(x, lambda, 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), @@ -1124,6 +1179,62 @@ test_that("dNmixture_BBP_s works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing values + # Uncompiled calculation + x <- c(1, 0, NA, 3, 0) + probX <- dNmixture_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + + # Manually calculate the correct answer + alpha <- prob * s + beta <- s - prob * s + + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dpois(N, lambda) * + prod(dBetaBinom_s(x, N, + alpha, shape2 = beta), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixture_BBP_s(x, lambda, prob, s, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Compilation and compiled calculations + CdNmixture_BBP_s <- compileNimble(dNmixture_BBP_s) + CprobX <- CdNmixture_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixture_BBP_s(x, lambda, 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), + 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_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixture_BBP_s(x, lambda, 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), @@ -1263,6 +1374,15 @@ test_that("dNmixture_BBP_oneObs works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing value + x <- as.numeric(NA) + probX <- dNmixture_BBP_oneObs(x, lambda, prob, s, Nmin, Nmax) + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixture_BBP_oneObs(x, lambda, prob, s, Nmin, Nmax) + expect_equal(CprobX, 1) + # Test imputing value for all NAs xNA <- NA mNA <- nimbleModel(nc, data = list(x = xNA), From d5563ec5b13aaab9c5abfaa89b5b19be072168b7 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 30 Aug 2024 11:12:07 -0400 Subject: [PATCH 07/11] Remove leftover code --- R/utils.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/utils.R b/R/utils.R index 920f96f..7eb2b5e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -191,8 +191,6 @@ dNmixture_BNB_steps <- nimbleFunction( numN <- NmaxC - NminC + 1 - 1 # remember... prods <- rep(0, numN) for (i in (NminC + 1):NmaxC){ - prods[i - NminC] <- (i + r - 1) * prod(i/(i - x)) / i - prodi <- i + r - 1 for (j in 1:length(x)){ xj <- ADbreak(x[j]) From ccaf0933738ca50133ad68a6a1d326848c82af0a Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 30 Aug 2024 11:44:30 -0400 Subject: [PATCH 08/11] dNmixture_BBP AD missing value support --- R/dNmixtureAD.R | 42 ++++-- tests/testthat/test-AD.R | 42 ++++++ tests/testthat/test-NmixtureADnoDerivs.R | 158 +++++++++++++++++++++++ 3 files changed, 234 insertions(+), 8 deletions(-) diff --git a/R/dNmixtureAD.R b/R/dNmixtureAD.R index 0ae5f09..aba812a 100644 --- a/R/dNmixtureAD.R +++ b/R/dNmixtureAD.R @@ -352,13 +352,26 @@ dNmixtureAD_BBP_v <- nimbleFunction( if (Nmin == -1 | Nmax == -1) { stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.") } - Nmin <- ADbreak(max( max(x), Nmin )) ## set Nmin to at least the largest x - logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, - dBetaBinom_v(x, Nmin, alpha, beta, len = len, log = TRUE)) + max_x <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + any_not_na <- TRUE + } + } + Nmin <- ADbreak(max( max_x, Nmin )) ## set Nmin to at least the largest x + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, + dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE), usingAD=TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs = list(run = list()) + }, buildDerivs = list(run = list(ignore=c("i", "xi"))) ) #' @rdname dNmixtureAD @@ -384,13 +397,26 @@ dNmixtureAD_BBP_s <- nimbleFunction( } #Clen <- 0L #Clen <- ADbreak(len) - Nmin <- ADbreak(max( max(x), Nmin )) ## set Nmin to at least the largest x - logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, - dBetaBinom_s(x, Nmin, alpha, beta, len = len, log = TRUE)) + max_x <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + any_not_na <- TRUE + } + } + Nmin <- ADbreak(max( max_x, Nmin )) ## set Nmin to at least the largest x + + logProb <- 0 + if(any_not_na){ + logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax, + dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE), usingAD=TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs = list(run=list()) + }, buildDerivs = list(run=list(ignore=c("i", "xi"))) ) #' @rdname dNmixtureAD diff --git a/tests/testthat/test-AD.R b/tests/testthat/test-AD.R index 6b7cb18..e042060 100644 --- a/tests/testthat/test-AD.R +++ b/tests/testthat/test-AD.R @@ -293,6 +293,27 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # Missing value handling + x <- c(7, 7, NA, 9, 10) + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, + lambda = lambda, + s = s), + buildDerivs = TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's')) + v1_case1 <- list(arg1 = c(prob, lambda, s)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ############################## #### dNmixture_BBNB_s case #### @@ -494,6 +515,27 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # Missing value handling + x <- c(7, 7, NA, 9, 10) + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, + lambda = lambda, + s = s), + buildDerivs=TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's')) + v1_case1 <- list(arg1 = c(prob, lambda, s)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ############################## #### dNmixture_BBNB_v case #### diff --git a/tests/testthat/test-NmixtureADnoDerivs.R b/tests/testthat/test-NmixtureADnoDerivs.R index c32a3da..ca2b425 100644 --- a/tests/testthat/test-NmixtureADnoDerivs.R +++ b/tests/testthat/test-NmixtureADnoDerivs.R @@ -1176,6 +1176,75 @@ test_that("dNmixtureAD_BBP_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) + Nmax <- 250 + probX <- dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) + + # Manually calculate the correct answer + alpha <- prob * s + beta <- s - prob * s + + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dpois(N, lambda) * + prod(dBetaBinom_v(x, N, shape1 = alpha, shape2 = beta), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Other Nmin / Nmax + Nmin <- 3 + for(Nmax in 3:6) { + dynProbX <- dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = Nmin, Nmax = Nmax, len) + dynCorrectProbX <- 0 + for (N in Nmin:Nmax) { + dynCorrectProbX <- dynCorrectProbX + dpois(N, lambda) * + prod(dBetaBinom_v(x, N, shape1 = alpha, shape2 = beta)) + } + expect_equal(dynProbX, dynCorrectProbX) + } + Nmin <- 0 + Nmax <- 250 + # Compilation and compiled calculations + CprobX <- CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_BBP_v(x, lambda, 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), + 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 <- dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixtureAD_BBP_v(x, lambda, 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), @@ -1347,6 +1416,77 @@ test_that("dNmixtureAD_BBP_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 <- dNmixtureAD_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + + # Manually calculate the correct answer + alpha <- prob * s + beta <- s - prob * s + + correctProbX <- 0 + for (N in Nmin:Nmax) { + correctProbX <- correctProbX + dpois(N, lambda) * + prod(dBetaBinom_s(x, N, + shape1 = alpha, shape2 = beta), na.rm=TRUE) + } + + expect_equal(probX, correctProbX) + + # Uncompiled log probability + lProbX <- dNmixtureAD_BBP_s(x, lambda, prob, s, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Other Nmin / Nmax + Nmin <- 3 + for(Nmax in 3:6) { + dynProbX <- dNmixtureAD_BBP_s(x, lambda, prob, s, Nmin = Nmin, Nmax = Nmax, len) + dynCorrectProbX <- 0 + for (N in Nmin:Nmax) { + dynCorrectProbX <- dynCorrectProbX + dpois(N, lambda) * + prod(dBetaBinom_s(x, N, shape1 = alpha, shape2 = beta)) + } + expect_equal(dynProbX, dynCorrectProbX) + } + Nmin <- 0 + Nmax <- 250 + + # Compilation and compiled calculations + CprobX <- CdNmixtureAD_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_BBP_s(x, lambda, 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), + 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 <- dNmixtureAD_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixtureAD_BBP_s(x, lambda, 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), @@ -1515,6 +1655,24 @@ test_that("dNmixtureAD_BBP_oneObs works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing value + # Uncompiled calculation + x <- as.numeric(NA) + lambda <- 8 + s <- 2 + prob <- c(0.5) + Nmin <- 0 + Nmax <- 250 + len <- 5 + + probX <- dNmixtureAD_BBP_oneObs(x, lambda, prob, s, Nmin, Nmax) + + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixtureAD_BBP_oneObs(x, lambda, prob, s, Nmin, Nmax) + expect_equal(CprobX, 1) + # Test imputing value for all NAs xNA <- NA mNA <- nimbleModel(nc, data = list(x = xNA), From 068cd1d83a73ccb0ea34201f4858d57e5c6f0486 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 30 Aug 2024 11:58:41 -0400 Subject: [PATCH 09/11] Code cleanup and fix a couple tests --- R/dNmixture.R | 18 ++++-------------- tests/testthat/test-Nmixture.R | 4 ++-- 2 files changed, 6 insertions(+), 16 deletions(-) diff --git a/R/dNmixture.R b/R/dNmixture.R index 6659aff..1680676 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -603,19 +603,14 @@ dNmixture_BBP_v <- nimbleFunction( stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.") } max_x <- 0 - for (i in 1:length(x)){ - if(!is.na(x[i])){ - if(x[i] > max_x) max_x <- x[i] - } - } - Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x - 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){ @@ -655,19 +650,14 @@ dNmixture_BBP_s <- nimbleFunction( } max_x <- 0 - for (i in 1:length(x)){ - if(!is.na(x[i])){ - if(x[i] > max_x) max_x <- x[i] - } - } - Nmin <- max( max_x, Nmin ) ## set Nmin to at least the largest x - 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){ diff --git a/tests/testthat/test-Nmixture.R b/tests/testthat/test-Nmixture.R index 4e467e2..6b51a71 100644 --- a/tests/testthat/test-Nmixture.R +++ b/tests/testthat/test-Nmixture.R @@ -1029,11 +1029,11 @@ test_that("dNmixture_BBP_v works", expect_equal(CMlProbX, lProbX) x <- as.numeric(rep(NA, 5)) - probX <- dNmixture_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + probX <- dNmixture_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) expect_equal(probX, 1) # Compilation and compiled calculations - CprobX <- CdNmixture_BBP_s(x, lambda, prob, s, Nmin, Nmax, len) + CprobX <- CdNmixture_BBP_v(x, lambda, prob, s, Nmin, Nmax, len) expect_equal(CprobX, 1) # Test imputing value for all NAs From 23dc177b3a0cd7708eb3141b45702e40c39d6815 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Tue, 3 Sep 2024 12:21:45 -0400 Subject: [PATCH 10/11] dNmixture BBNB (non-AD) missing value support --- R/dNmixture.R | 38 ++++++++-- R/utils.R | 15 ++-- tests/testthat/test-Nmixture.R | 125 +++++++++++++++++++++++++++++++++ 3 files changed, 168 insertions(+), 10 deletions(-) diff --git a/R/dNmixture.R b/R/dNmixture.R index 1680676..9df0069 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -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()) @@ -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()) diff --git a/R/utils.R b/R/utils.R index 7eb2b5e..3ce140b 100644 --- a/R/utils.R +++ b/R/utils.R @@ -303,9 +303,16 @@ 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 @@ -313,5 +320,5 @@ dNmixture_BBNB_steps <- nimbleFunction( return(logProb) returnType(double()) }, - buildDerivs = list(run = list(ignore = c("i"))) + buildDerivs = list(run = list(ignore = c("i","j","xj"))) ) diff --git a/tests/testthat/test-Nmixture.R b/tests/testthat/test-Nmixture.R index 6b51a71..e86c4c7 100644 --- a/tests/testthat/test-Nmixture.R +++ b/tests/testthat/test-Nmixture.R @@ -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), @@ -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), @@ -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), From f34b98e805a45010fc34f7d0b22d172cdbe62a0c Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Tue, 10 Sep 2024 14:26:08 -0400 Subject: [PATCH 11/11] BBNB AD missing value support --- R/dNmixtureAD.R | 43 +++++-- tests/testthat/test-AD.R | 62 +++++++++ tests/testthat/test-NmixtureADnoDerivs.R | 157 +++++++++++++++++++++++ 3 files changed, 254 insertions(+), 8 deletions(-) diff --git a/R/dNmixtureAD.R b/R/dNmixtureAD.R index aba812a..beb70eb 100644 --- a/R/dNmixtureAD.R +++ b/R/dNmixtureAD.R @@ -462,13 +462,27 @@ dNmixtureAD_BBNB_v <- nimbleFunction( if (Nmin == -1 | Nmax == -1) { stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.") } - Nmin <- ADbreak(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 = len, log = TRUE)) + + max_x <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + any_not_na <- TRUE + } + } + Nmin <- ADbreak(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), usingAD=TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs=list(run=list()) + }, buildDerivs=list(run=list(ignore=c("i", "xi"))) ) #' @rdname dNmixtureAD @@ -500,13 +514,26 @@ dNmixtureAD_BBNB_s <- nimbleFunction( } # Clen <- 0L # Clen <- ADbreak(len) - Nmin <- ADbreak(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 = len, log = TRUE)) + max_x <- 0 + any_not_na <- FALSE + for (i in 1:length(x)){ + xi <- ADbreak(x[i]) + if(!is.na(xi)){ + if(x[i] > max_x) max_x <- x[i] + any_not_na <- TRUE + } + } + Nmin <- ADbreak(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), usingAD=TRUE) + } if (log) return(logProb) else return(exp(logProb)) returnType(double()) - }, buildDerivs=list(run=list()) + }, buildDerivs=list(run=list(ignore=c("i", "xi"))) ) #' @rdname dNmixtureAD diff --git a/tests/testthat/test-AD.R b/tests/testthat/test-AD.R index e042060..890047e 100644 --- a/tests/testthat/test-AD.R +++ b/tests/testthat/test-AD.R @@ -354,6 +354,28 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + # Missing value + x <- c(7, 7, NA, 9, 10) + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, + lambda = lambda, + theta = theta, s = s), + buildDerivs=TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's')) + v1_case1 <- list(arg1 = c(prob, lambda, theta, s)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + + ########################## #### dNmixture_v case #### @@ -580,6 +602,26 @@ test_that ("dNmixture works with AD", { v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + x <- c(7, 7, NA, 9, 10) + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, + lambda = lambda, + theta = theta, s = s), + buildDerivs=TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's')) + v1_case1 <- list(arg1 = c(prob, lambda, theta, s)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + ########################## #### dNmixture_BNB_oneObs case #### @@ -692,6 +734,26 @@ test_that ("dNmixture works with AD", { model_calculate_test, nodesList_case1, v1_case1, v2_case1, 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + + x <- as.numeric(NA) + Rmodel <- nimbleModel(nc, data = list(x = x), + inits = list(prob = prob, theta = theta, s=s, + lambda = lambda), + buildDerivs=TRUE) + Rmodel$calculate() + + Cmodel <- compileNimble(Rmodel) + Cmodel$calculate() + + nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's')) + v1_case1 <- list(arg1 = c(prob, lambda, theta, s2)) # taping values for prob and lambda + v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda + + res <- model_calculate_test_case(Rmodel, Cmodel, + model_calculate_test, nodesList_case1, + v1_case1, v2_case1, + 0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14)) + }) diff --git a/tests/testthat/test-NmixtureADnoDerivs.R b/tests/testthat/test-NmixtureADnoDerivs.R index ca2b425..137edcb 100644 --- a/tests/testthat/test-NmixtureADnoDerivs.R +++ b/tests/testthat/test-NmixtureADnoDerivs.R @@ -1848,6 +1848,84 @@ test_that("dNmixtureAD_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 <- dNmixtureAD_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 <- dNmixtureAD_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + + # Other Nmin / Nmax + Nmin <- 3 + for(Nmax in 3:6) { + dynProbX <- dNmixtureAD_BBNB_v(x, lambda, theta, prob, s, Nmin = Nmin, Nmax = Nmax, len) + dynCorrectProbX <- 0 + for (N in Nmin:Nmax) { + dynCorrectProbX <- dynCorrectProbX + dnbinom(N, size = r, prob = pNB) * + prod(dBetaBinom_v(x, N, shape1 = alpha, shape2 = beta)) + } + expect_equal(dynProbX, dynCorrectProbX) + } + Nmin <- 0 + Nmax <- 250 + + CprobX <- CdNmixtureAD_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE) + expect_equal(ClProbX, lProbX) + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixtureAD_BBNB_v(lambda = lambda, prob = prob[1:5], + theta = theta, s = s, + Nmin = Nmin, Nmax = Nmax, len = len) + + }) + + 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) + + # All missing + x <- as.numeric(rep(NA, 5)) + probX <- dNmixtureAD_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + CprobX <- CdNmixtureAD_BBNB_v(x, lambda, theta, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -2017,6 +2095,76 @@ test_that("dNmixtureAD_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 <- dNmixtureAD_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 <- dNmixtureAD_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len, log = TRUE) + lCorrectProbX <- log(correctProbX) + expect_equal(lProbX, lCorrectProbX) + # Other Nmin / Nmax + Nmin <- 3 + for(Nmax in 3:6) { + dynProbX <- dNmixtureAD_BBNB_s(x, lambda, theta, prob, s, Nmin = Nmin, Nmax = Nmax, len) + dynCorrectProbX <- 0 + for (N in Nmin:Nmax) { + dynCorrectProbX <- dynCorrectProbX + dnbinom(N, size = r, prob = pNB) * + prod(dBetaBinom_s(x, N, shape1 = alpha, shape2 = beta)) + } + expect_equal(dynProbX, dynCorrectProbX) + } + Nmin <- 0 + Nmax <- 250 + + CprobX <- CdNmixtureAD_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + + ClProbX <- CdNmixtureAD_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) + + # All missing + x <- as.numeric(rep(NA,5)) + probX <- dNmixtureAD_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + CprobX <- CdNmixtureAD_BBNB_s(x, lambda, theta, prob, s, Nmin, Nmax, len) + expect_equal(CprobX, probX) + # Test imputing value for all NAs xNA <- c(NA, NA, NA, NA, NA) mNA <- nimbleModel(nc, data = list(x = xNA), @@ -2186,6 +2334,15 @@ test_that("dNmixtureAD_BBNB_oneObs works", CMlProbX <- cm$getLogProb("x") expect_equal(CMlProbX, lProbX) + # Missing value + x <- as.numeric(NA) + probX <- dNmixtureAD_BBNB_oneObs(x, lambda, theta, prob, s, Nmin, Nmax) + expect_equal(probX, 1) + + CprobX <- CdNmixtureAD_BBNB_oneObs(x, lambda, theta, prob, s, Nmin, Nmax) + expect_equal(CprobX, probX) + + # Test imputing value for all NAs xNA <- NA mNA <- nimbleModel(nc, data = list(x = xNA),