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/R/dNmixture.R b/R/dNmixture.R index 4edaa23..9df0069 100644 --- a/R/dNmixture.R +++ b/R/dNmixture.R @@ -241,13 +241,50 @@ 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))) - 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))) + 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)){ + 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()) @@ -273,13 +310,50 @@ 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))) - 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))) + 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)){ + 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()) @@ -357,13 +431,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()) @@ -395,13 +508,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()) @@ -450,9 +602,21 @@ 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 + 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_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()) @@ -484,9 +648,22 @@ 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 + 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_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()) @@ -537,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()) @@ -575,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/dNmixtureAD.R b/R/dNmixtureAD.R index 76e13d8..beb70eb 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 @@ -163,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 @@ -190,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 @@ -255,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 @@ -287,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 @@ -339,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 @@ -377,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/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/R/utils.R b/R/utils.R index ef511ee..3ce140b 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 ##### @@ -183,8 +190,16 @@ dNmixture_BNB_steps <- nimbleFunction( numN <- 0L 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 + for (i in (NminC + 1):NmaxC){ + 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 @@ -192,7 +207,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"))) ) @@ -230,8 +245,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 @@ -239,7 +262,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"))) ) @@ -280,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 @@ -290,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-AD.R b/tests/testthat/test-AD.R index 3089fc6..890047e 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", { @@ -124,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 #### @@ -162,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 #### @@ -200,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 #### @@ -240,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 #### @@ -278,6 +414,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 #### @@ -320,6 +474,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 #### @@ -362,6 +537,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 #### @@ -406,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 #### @@ -518,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-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), diff --git a/tests/testthat/test-Nmixture.R b/tests/testthat/test-Nmixture.R index 32f1e81..e86c4c7 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), @@ -372,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), @@ -506,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), @@ -639,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), @@ -782,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_v(x, lambda, prob, s, Nmin, Nmax, len) + expect_equal(probX, 1) + + # Compilation and compiled calculations + CprobX <- CdNmixture_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), @@ -925,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), @@ -1064,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), @@ -1209,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), @@ -1351,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), @@ -1494,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), diff --git a/tests/testthat/test-NmixtureADnoDerivs.R b/tests/testthat/test-NmixtureADnoDerivs.R index 67fe741..137edcb 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), @@ -388,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), @@ -543,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), @@ -692,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), @@ -805,43 +1122,104 @@ test_that("dNmixtureAD_BBP_v works", return(dNmixtureAD_BBP_v(x,lambda,prob,s,Nmin,Nmax,len,log)) returnType(double()) } - ) + ) + # Compilation and compiled calculations + CdNmixtureAD_BBP_v <- compileNimble(call_dNmixtureAD_BBP_v) + 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) + + # Dynamic Nmin / Nmax isn't allowed + expect_error({ + dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax = -1, len) + }) + expect_error({ + dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax, len) + }) + expect_error({ + dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax = -1, len) + }) + expect_error({ + CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax = -1, len) + }) + expect_error({ + CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax, len) + }) + expect_error({ + CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax = -1, len) + }) + + + # Use in Nimble model + nc <- nimbleCode({ + x[1:5] ~ dNmixtureAD_BBP_v(lambda = lambda, prob = prob[1:5], + s = s, + Nmin = Nmin, Nmax = Nmax, len = len) + }) + + 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) + + # 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 - CdNmixtureAD_BBP_v <- compileNimble(call_dNmixtureAD_BBP_v) 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) - # Dynamic Nmin / Nmax isn't allowed - expect_error({ - dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax = -1, len) - }) - expect_error({ - dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax, len) - }) - expect_error({ - dNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax = -1, len) - }) - expect_error({ - CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax = -1, len) - }) - expect_error({ - CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax, len) - }) - expect_error({ - CdNmixtureAD_BBP_v(x, lambda, prob, s, Nmin, Nmax = -1, len) - }) - - - # Use in Nimble model - nc <- nimbleCode({ - x[1:5] ~ dNmixtureAD_BBP_v(lambda = lambda, prob = prob[1:5], - s = s, - Nmin = Nmin, Nmax = Nmax, len = len) - }) - m <- nimbleModel(code = nc, data = list(x = x), inits = list(lambda = lambda, @@ -859,6 +1237,14 @@ test_that("dNmixtureAD_BBP_v works", 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), @@ -1030,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), @@ -1198,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), @@ -1373,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), @@ -1542,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), @@ -1711,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), 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)