Skip to content

Commit

Permalink
Merge pull request #29 from kenkellner/handle_na
Browse files Browse the repository at this point in the history 
Handle missing values in `dOcc`* and `dNmixture`* distributions
  • Loading branch information
@dochvam
dochvam authored Oct 1, 2024
2 parents 38e9e88 + f34b98e commit 5b6c41d
Show file tree
Hide file tree
Showing 10 changed files with 1,931 additions and 132 deletions.
26 changes: 16 additions & 10 deletions R/dBetaBinom.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
283 changes: 243 additions & 40 deletions R/dNmixture.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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())
Expand All @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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())
Expand Down
Loading

0 comments on commit 5b6c41d

Please sign in to comment.