-
Notifications
You must be signed in to change notification settings - Fork 2
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Instead of fixing bottleneck size for CompartmentType, draw from distribution #126
Comments
I think this can be accomplished by modifying this line: Line 244 in ae05f46
If we let the class mean.size <- comp$get.type()$get.bottleneck.size()
bottleneck.size <- rpois(n=1, lambda=mean.size) For example, this lets us assign a mean bottleneck size of 1.2 to a CompartmentType, and then the realized bottleneck size for a given event has some probability of being 1, 2 or more, conditional on being more than zero. |
We want to use a zero-truncated Poisson or negative binomial distribution: |
We should use a negative binomial distribution because the number of transmitted-founders is probably not Poisson-distributed: |
https://rdrr.io/rforge/countreg/man/ztnbinom.html ## ztnbinom: Zero-truncated negative binomial
dztnbinom <- function(x, mu, theta, size, log = FALSE) {
if(!missing(theta) & !missing(size)) stop("only 'theta' or 'size' may be specified")
if(!missing(size)) theta <- size
rval <- dnbinom(x, mu = mu, size = theta, log = TRUE) - pnbinom(0, mu = mu, size = theta, lower.tail = FALSE, log.p = TRUE)
rval[x < 1] <- -Inf
rval[mu <= 0] <- 0
if(log) rval else exp(rval)
}
pztnbinom <- function(q, mu, theta, size, lower.tail = TRUE, log.p = FALSE) {
if(!missing(theta) & !missing(size)) stop("only 'theta' or 'size' may be specified")
if(!missing(size)) theta <- size
rval <- log(pnbinom(q, mu = mu, size = theta, lower.tail = lower.tail, log.p = FALSE) - dnbinom(0, mu = mu, size = theta)) -
pnbinom(0, mu = mu, size = theta, lower.tail = FALSE, log.p = TRUE)
rval[q < 1] <- if(lower.tail) -Inf else 0
if(log.p) rval else exp(rval)
}
qztnbinom <- function(p, mu, theta, size, lower.tail = TRUE, log.p = FALSE) {
if(!missing(theta) & !missing(size)) stop("only 'theta' or 'size' may be specified")
if(!missing(size)) theta <- size
p_orig <- p
p <- if(log.p) p else log(p)
p <- p + pnbinom(0, mu = mu, size = theta, lower.tail = FALSE, log.p = TRUE)
p <- exp(p) + dnbinom(0, mu = mu, size = theta)
rval <- qnbinom(p, mu = mu, size = theta, lower.tail = lower.tail, log.p = FALSE)
if(lower.tail) rval[p_orig < dztnbinom(1, mu = mu, theta = theta, log = log.p)] <- 1
rval
}
rztnbinom <- function(n, mu, theta, size) {
if(!missing(theta) & !missing(size)) stop("only 'theta' or 'size' may be specified")
if(!missing(size)) theta <- size
qztnbinom(runif(n), mu = mu, theta = theta)
} |
This needs testing before we can close this issue |
So each Compartment of that Type is assigned a random deviate from the distribution, e.g., Poisson with Type-specific mean \lambda.
The text was updated successfully, but these errors were encountered: