diff --git a/DESCRIPTION b/DESCRIPTION index eb1b244c..4975d5b1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.19 +Version: 0.2.20 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), diff --git a/R/cpp11.R b/R/cpp11.R index 9b37f389..15a17c9c 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -40,8 +40,12 @@ monty_rng_binomial <- function(ptr, n, r_size, r_prob, n_threads, is_float) { .Call(`_monty_monty_rng_binomial`, ptr, n, r_size, r_prob, n_threads, is_float) } -monty_rng_nbinomial <- function(ptr, n, r_size, r_prob, n_threads, is_float) { - .Call(`_monty_monty_rng_nbinomial`, ptr, n, r_size, r_prob, n_threads, is_float) +monty_rng_negative_binomial_prob <- function(ptr, n, r_size, r_prob, n_threads, is_float) { + .Call(`_monty_monty_rng_negative_binomial_prob`, ptr, n, r_size, r_prob, n_threads, is_float) +} + +monty_rng_negative_binomial_mu <- function(ptr, n, r_size, r_mu, n_threads, is_float) { + .Call(`_monty_monty_rng_negative_binomial_mu`, ptr, n, r_size, r_mu, n_threads, is_float) } monty_rng_hypergeometric <- function(ptr, n, r_n1, r_n2, r_k, n_threads, is_float) { diff --git a/R/dsl-distributions.R b/R/dsl-distributions.R index d8af6aae..6b9ee832 100644 --- a/R/dsl-distributions.R +++ b/R/dsl-distributions.R @@ -103,7 +103,7 @@ distr_gamma_rate <- distribution( distr_gamma_scale <- distribution( name = "Gamma", - variant = "rate", + variant = "scale", density = function(x, shape, scale) { dgamma(x, shape, scale = scale, log = TRUE) }, @@ -125,6 +125,36 @@ distr_hypergeometric <- distribution( mean = quote(k * n1 / (n1 + n2))), cpp = list(density = "hypergeometric", sample = "hypergeometric")) +distr_negative_binomial_prob <- distribution( + name = "NegativeBinomial", + variant = "prob", + density = function(x, size, prob) { + dnbinom(x, size, prob = prob, log = TRUE) + }, + domain = c(0, Inf), + sample = function(rng, size, prob) rng$negative_binomial_prob(1, size, prob), + expr = list( + density = quote(lgamma(x + size) - lgamma(size) - lgamma(x + 1) + + x * log(1 - prob) + size * log(prob)), + mean = quote(size * (1 - prob) / prob)), + cpp = list(density = "negative_binomial_prob", + sample = "negative_binomial_prob")) + +distr_negative_binomial_mu <- distribution( + name = "NegativeBinomial", + variant = "mu", + density = function(x, size, mu) { + dnbinom(x, size, mu = mu, log = TRUE) + }, + domain = c(0, Inf), + sample = function(rng, size, mu) rng$negative_binomial_mu(1, size, mu), + expr = list( + density = quote(lgamma(x + size) - lgamma(size) - lgamma(x + 1) + + size * log(size) + x * log(mu) - + (size + x) * log(size + mu)), + mean = quote(mu)), + cpp = list(density = "negative_binomial_mu", sample = "negative_binomial_mu")) + distr_normal <- distribution( name = "Normal", density = function(x, mean, sd) dnorm(x, mean, sd, log = TRUE), @@ -166,6 +196,8 @@ dsl_distributions <- local({ distr_gamma_rate, distr_gamma_scale, distr_hypergeometric, + distr_negative_binomial_prob, + distr_negative_binomial_mu, distr_normal, distr_poisson, distr_uniform) diff --git a/R/rng.R b/R/rng.R index 35c98919..27cf60d7 100644 --- a/R/rng.R +++ b/R/rng.R @@ -96,7 +96,10 @@ ##' rng$binomial(5, 10, 0.3) ##' ##' # Negative binomially distributed random numbers with size and prob -##' rng$nbinomial(5, 10, 0.3) +##' rng$negative_binomial_prob(5, 10, 0.3) +##' +##' # Negative binomially distributed random numbers with size and mean mu +##' rng$negative_binomial_mu(5, 10, 25) ##' ##' # Hypergeometric distributed random numbers with parameters n1, n2 and k ##' rng$hypergeometric(5, 6, 10, 4) @@ -283,9 +286,25 @@ monty_rng <- R6::R6Class( ##' (between 0 and 1, length 1 or n) ##' ##' @param n_threads Number of threads to use; see Details - nbinomial = function(n, size, prob, n_threads = 1L) { - monty_rng_nbinomial(private$ptr, n, size, prob, n_threads, - private$float) + negative_binomial_prob = function(n, size, prob, n_threads = 1L) { + monty_rng_negative_binomial_prob(private$ptr, n, size, prob, n_threads, + private$float) + }, + + ##' @description Generate `n` numbers from a negative binomial distribution + ##' + ##' @param n Number of samples to draw (per stream) + ##' + ##' @param size The target number of successful trials + ##' (zero or more, length 1 or n) + ##' + ##' @param mu The mean + ##' (zero or more, length 1 or n) + ##' + ##' @param n_threads Number of threads to use; see Details + negative_binomial_mu = function(n, size, mu, n_threads = 1L) { + monty_rng_negative_binomial_mu(private$ptr, n, size, mu, n_threads, + private$float) }, ##' @description Generate `n` numbers from a hypergeometric distribution diff --git a/inst/include/monty/random/nbinomial.hpp b/inst/include/monty/random/nbinomial.hpp deleted file mode 100644 index 004623fb..00000000 --- a/inst/include/monty/random/nbinomial.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#pragma once - -#include - -#include "monty/random/gamma.hpp" -#include "monty/random/poisson.hpp" -#include "monty/random/generator.hpp" - -namespace monty { -namespace random { - -namespace { - -template -void nbinomial_validate(real_type size, real_type prob) { - if(!R_FINITE(size) || !R_FINITE(prob) || size <= 0 || prob <= 0 || prob > 1) { - char buffer[256]; - snprintf(buffer, 256, - "Invalid call to nbinomial with size = %g, prob = %g", - size, prob); - monty::utils::fatal_error(buffer); - } -} - -} - -template -real_type nbinomial(rng_state_type& rng_state, real_type size, real_type prob) { -#ifdef __CUDA_ARCH__ - static_assert("nbinomial() not implemented for GPU targets"); -#endif - nbinomial_validate(size, prob); - - if (rng_state.deterministic) { - return (1 - prob) * size / prob; - } - return (prob == 1) ? 0 : poisson(rng_state, gamma_scale(rng_state, size, (1 - prob) / prob)); -} - -} -} diff --git a/inst/include/monty/random/negative_binomial.hpp b/inst/include/monty/random/negative_binomial.hpp new file mode 100644 index 00000000..2c3cb61d --- /dev/null +++ b/inst/include/monty/random/negative_binomial.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include + +#include "monty/random/gamma.hpp" +#include "monty/random/poisson.hpp" +#include "monty/random/generator.hpp" + +namespace monty { +namespace random { + +namespace { + +template +void negative_binomial_validate(real_type size, real_type prob) { + if(!R_FINITE(size) || !R_FINITE(prob) || size <= 0 || prob <= 0 || prob > 1) { + char buffer[256]; + snprintf(buffer, 256, + "Invalid call to negative_binomial with size = %g, prob = %g", + size, prob); + monty::utils::fatal_error(buffer); + } +} + +} + +template +real_type negative_binomial_prob(rng_state_type& rng_state, real_type size, real_type prob) { +#ifdef __CUDA_ARCH__ + static_assert("negative_binomial_prob() not implemented for GPU targets"); +#endif + negative_binomial_validate(size, prob); + + if (rng_state.deterministic) { + return (1 - prob) * size / prob; + } + return (prob == 1) ? 0 : poisson(rng_state, gamma_scale(rng_state, size, (1 - prob) / prob)); +} + +template +real_type negative_binomial_mu(rng_state_type& rng_state, real_type size, real_type mu) { +#ifdef __CUDA_ARCH__ + static_assert("negative_binomial_mu() not implemented for GPU targets"); +#endif + const auto prob = size / (size + mu); + negative_binomial_validate(size, prob); + + return negative_binomial_prob(rng_state, size, prob); +} + +} +} diff --git a/inst/include/monty/random/random.hpp b/inst/include/monty/random/random.hpp index 508a52d8..93df6424 100644 --- a/inst/include/monty/random/random.hpp +++ b/inst/include/monty/random/random.hpp @@ -10,7 +10,7 @@ #include "monty/random/gamma.hpp" #include "monty/random/hypergeometric.hpp" #include "monty/random/multinomial.hpp" -#include "monty/random/nbinomial.hpp" +#include "monty/random/negative_binomial.hpp" #include "monty/random/normal.hpp" #include "monty/random/poisson.hpp" #include "monty/random/uniform.hpp" diff --git a/man/monty_rng.Rd b/man/monty_rng.Rd index 6962e231..4642f365 100644 --- a/man/monty_rng.Rd +++ b/man/monty_rng.Rd @@ -97,7 +97,10 @@ rng$normal(5, 4, 2) rng$binomial(5, 10, 0.3) # Negative binomially distributed random numbers with size and prob -rng$nbinomial(5, 10, 0.3) +rng$negative_binomial_prob(5, 10, 0.3) + +# Negative binomially distributed random numbers with size and mean mu +rng$negative_binomial_mu(5, 10, 25) # Hypergeometric distributed random numbers with parameters n1, n2 and k rng$hypergeometric(5, 6, 10, 4) @@ -140,7 +143,8 @@ rng$multinomial(5, 10, c(0.1, 0.3, 0.5, 0.1)) \item \href{#method-monty_rng-uniform}{\code{monty_rng$uniform()}} \item \href{#method-monty_rng-normal}{\code{monty_rng$normal()}} \item \href{#method-monty_rng-binomial}{\code{monty_rng$binomial()}} -\item \href{#method-monty_rng-nbinomial}{\code{monty_rng$nbinomial()}} +\item \href{#method-monty_rng-negative_binomial_prob}{\code{monty_rng$negative_binomial_prob()}} +\item \href{#method-monty_rng-negative_binomial_mu}{\code{monty_rng$negative_binomial_mu()}} \item \href{#method-monty_rng-hypergeometric}{\code{monty_rng$hypergeometric()}} \item \href{#method-monty_rng-gamma_scale}{\code{monty_rng$gamma_scale()}} \item \href{#method-monty_rng-gamma_rate}{\code{monty_rng$gamma_rate()}} @@ -341,12 +345,12 @@ Generate \code{n} numbers from a binomial distribution } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-nbinomial}{}}} -\subsection{Method \code{nbinomial()}}{ +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-monty_rng-negative_binomial_prob}{}}} +\subsection{Method \code{negative_binomial_prob()}}{ Generate \code{n} numbers from a negative binomial distribution \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$nbinomial(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{monty_rng$negative_binomial_prob(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -360,6 +364,31 @@ Generate \code{n} numbers from a negative binomial distribution \item{\code{prob}}{The probability of success on each trial (between 0 and 1, length 1 or n)} +\item{\code{n_threads}}{Number of threads to use; see Details} +} +\if{html}{\out{}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-monty_rng-negative_binomial_mu}{}}} +\subsection{Method \code{negative_binomial_mu()}}{ +Generate \code{n} numbers from a negative binomial distribution +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{monty_rng$negative_binomial_mu(n, size, mu, n_threads = 1L)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{n}}{Number of samples to draw (per stream)} + +\item{\code{size}}{The target number of successful trials +(zero or more, length 1 or n)} + +\item{\code{mu}}{The mean +(zero or more, length 1 or n)} + \item{\code{n_threads}}{Number of threads to use; see Details} } \if{html}{\out{
}} diff --git a/src/cpp11.cpp b/src/cpp11.cpp index d2dc2697..b5943761 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -78,10 +78,17 @@ extern "C" SEXP _monty_monty_rng_binomial(SEXP ptr, SEXP n, SEXP r_size, SEXP r_ END_CPP11 } // random.cpp -cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads, bool is_float); -extern "C" SEXP _monty_monty_rng_nbinomial(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads, SEXP is_float) { +cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads, bool is_float); +extern "C" SEXP _monty_monty_rng_negative_binomial_prob(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads, SEXP is_float) { BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_nbinomial(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(is_float))); + return cpp11::as_sexp(monty_rng_negative_binomial_prob(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(is_float))); + END_CPP11 +} +// random.cpp +cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_mu, int n_threads, bool is_float); +extern "C" SEXP _monty_monty_rng_negative_binomial_mu(SEXP ptr, SEXP n, SEXP r_size, SEXP r_mu, SEXP n_threads, SEXP is_float) { + BEGIN_CPP11 + return cpp11::as_sexp(monty_rng_negative_binomial_mu(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_mu), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(is_float))); END_CPP11 } // random.cpp @@ -208,7 +215,8 @@ static const R_CallMethodDef CallEntries[] = { {"_monty_monty_rng_jump", (DL_FUNC) &_monty_monty_rng_jump, 2}, {"_monty_monty_rng_long_jump", (DL_FUNC) &_monty_monty_rng_long_jump, 2}, {"_monty_monty_rng_multinomial", (DL_FUNC) &_monty_monty_rng_multinomial, 6}, - {"_monty_monty_rng_nbinomial", (DL_FUNC) &_monty_monty_rng_nbinomial, 6}, + {"_monty_monty_rng_negative_binomial_mu", (DL_FUNC) &_monty_monty_rng_negative_binomial_mu, 6}, + {"_monty_monty_rng_negative_binomial_prob", (DL_FUNC) &_monty_monty_rng_negative_binomial_prob, 6}, {"_monty_monty_rng_normal", (DL_FUNC) &_monty_monty_rng_normal, 7}, {"_monty_monty_rng_pointer_init", (DL_FUNC) &_monty_monty_rng_pointer_init, 4}, {"_monty_monty_rng_pointer_sync", (DL_FUNC) &_monty_monty_rng_pointer_sync, 2}, diff --git a/src/random.cpp b/src/random.cpp index cc16b0e4..88e95cf8 100644 --- a/src/random.cpp +++ b/src/random.cpp @@ -336,7 +336,7 @@ cpp11::sexp monty_rng_binomial(SEXP ptr, int n, template -cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, +cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads) { T *rng = cpp11::as_cpp>(ptr).get(); @@ -363,7 +363,7 @@ cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, for (size_t j = 0; j < (size_t)n; ++j) { auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; auto prob_ij = prob_vary.draw ? prob_i[j] : prob_i[0]; - y_i[j] = monty::random::nbinomial(state, size_ij, prob_ij); + y_i[j] = monty::random::negative_binomial_prob(state, size_ij, prob_ij); } } catch (std::exception const& e) { errors.capture(e, i); @@ -375,6 +375,48 @@ cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, return sexp_matrix(ret, n, n_streams); } + +template +cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, + cpp11::doubles r_size, cpp11::doubles r_mu, + int n_threads) { + T *rng = cpp11::as_cpp>(ptr).get(); + const int n_streams = rng->size(); + cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); + double * y = REAL(ret); + + const double * size = REAL(r_size); + const double * mu = REAL(r_mu); + auto size_vary = check_input_type(r_size, n, n_streams, "size"); + auto mu_vary = check_input_type(r_mu, n, n_streams, "mu"); + + monty::utils::openmp_errors errors(n_streams); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) num_threads(n_threads) +#endif + for (int i = 0; i < n_streams; ++i) { + try { + auto &state = rng->state(i); + auto y_i = y + n * i; + auto size_i = size_vary.generator ? size + size_vary.offset * i : size; + auto mu_i = mu_vary.generator ? mu + mu_vary.offset * i : mu; + for (size_t j = 0; j < (size_t)n; ++j) { + auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; + auto mu_ij = mu_vary.draw ? mu_i[j] : mu_i[0]; + y_i[j] = monty::random::negative_binomial_mu(state, size_ij, mu_ij); + } + } catch (std::exception const& e) { + errors.capture(e, i); + } + } + + errors.report("generators", 4, true); + + return sexp_matrix(ret, n, n_streams); +} + + template cpp11::sexp monty_rng_poisson(SEXP ptr, int n, cpp11::doubles r_lambda, int n_threads) { @@ -812,12 +854,21 @@ cpp11::sexp monty_rng_binomial(SEXP ptr, int n, } [[cpp11::register]] -cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, +cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads, bool is_float) { return is_float ? - monty_rng_nbinomial(ptr, n, r_size, r_prob, n_threads) : - monty_rng_nbinomial(ptr, n, r_size, r_prob, n_threads); + monty_rng_negative_binomial_prob(ptr, n, r_size, r_prob, n_threads) : + monty_rng_negative_binomial_prob(ptr, n, r_size, r_prob, n_threads); +} + +[[cpp11::register]] +cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, + cpp11::doubles r_size, cpp11::doubles r_mu, + int n_threads, bool is_float) { + return is_float ? + monty_rng_negative_binomial_mu(ptr, n, r_size, r_mu, n_threads) : + monty_rng_negative_binomial_mu(ptr, n, r_size, r_mu, n_threads); } [[cpp11::register]] diff --git a/tests/testthat/test-rng.R b/tests/testthat/test-rng.R index 41de599b..2b022aa4 100644 --- a/tests/testthat/test-rng.R +++ b/tests/testthat/test-rng.R @@ -1357,13 +1357,24 @@ test_that("can generate negative binomial numbers", { m <- 1000000 n <- 958 p <- 0.004145 - yf <- monty_rng$new(1)$nbinomial(m, n, p) + yf <- monty_rng$new(1)$negative_binomial_prob(m, n, p) expect_equal(mean(yf), (1 - p) * n / p, tolerance = 1e-3) expect_equal(var(yf), ((1 - p) * n) / p^2, tolerance = 1e-2) }) +test_that("negative_binomial_mu follows from negative_binomial_prob", { + rng1 <- monty_rng$new(seed = 1L) + rng2 <- monty_rng$new(seed = 1L) + prob <- 0.3 + size <- 20 + expect_identical( + rng1$negative_binomial_mu(100, size, size * (1 - prob) / prob), + rng2$negative_binomial_prob(100, size, prob)) +}) + + test_that("deterministic negative binomial returns mean", { m <- 100 p <- as.numeric(sample(10, m, replace = TRUE)) / 10 @@ -1372,22 +1383,23 @@ test_that("deterministic negative binomial returns mean", { rng_f <- monty_rng$new(1, real_type = "float", deterministic = TRUE) rng_d <- monty_rng$new(1, real_type = "double", deterministic = TRUE) - expect_equal(rng_f$nbinomial(m, n, p), (1 - p) * n / p, tolerance = 1e-6) - expect_equal(rng_d$nbinomial(m, n, p), (1 - p) * n / p) + expect_equal(rng_f$negative_binomial_prob(m, n, p), (1 - p) * n / p, + tolerance = 1e-6) + expect_equal(rng_d$negative_binomial_prob(m, n, p), (1 - p) * n / p) }) test_that("negative binomial prevents bad inputs", { - expect_error(monty_rng$new(1)$nbinomial(1, 10, 0), - "Invalid call to nbinomial with size = 10, prob = 0") - expect_error(monty_rng$new(1)$nbinomial(1, 0, 0.5), - "Invalid call to nbinomial with size = 0, prob = 0.5") - expect_error(monty_rng$new(1)$nbinomial(1, 10, 1.5), - "Invalid call to nbinomial with size = 10, prob = 1.5") - expect_error(monty_rng$new(1)$nbinomial(1, 10, Inf), - "Invalid call to nbinomial with size = 10, prob = inf") - expect_error(monty_rng$new(1)$nbinomial(1, Inf, 0.4), - "Invalid call to nbinomial with size = inf, prob = 0.4") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 10, 0), + "Invalid call to negative_binomial with size = 10, prob = 0") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 0, 0.5), + "Invalid call to negative_binomial with size = 0, prob = 0.5") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 10, 1.5), + "Invalid call to negative_binomial with size = 10, prob = 1.5") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 10, Inf), + "Invalid call to negative_binomial with size = 10, prob = inf") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, Inf, 0.4), + "Invalid call to negative_binomial with size = inf, prob = 0.4") })