diff --git a/DESCRIPTION b/DESCRIPTION index 51081a76..332cc5f7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: extras Title: Helper Functions for Bayesian Analyses -Version: 0.7.3.9002 +Version: 0.8.0 Authors@R: c( person("Nicole", "Hill", , "nicole@poissonconsulting.ca", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7623-2153")), @@ -32,6 +32,7 @@ Suggests: ggplot2, hms, knitr, + memoise, rlang, rmarkdown, scales, @@ -41,10 +42,11 @@ Suggests: tidyr, viridis, withr +VignetteBuilder: + knitr +Config/Needs/website: poissonconsulting/poissontemplate Config/testthat/edition: 3 Encoding: UTF-8 Language: en-US Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2 -VignetteBuilder: knitr -Config/Needs/website: poissonconsulting/poissontemplate +RoxygenNote: 7.3.2.9000 diff --git a/NEWS.md b/NEWS.md index 5fb4d984..992fd516 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,11 @@ +# extras 0.8.0 + +- Added a scalar case to `log_lik_beta_binom()` to improve speed for scalar inputs. +- Add memoization (if memoize package is installed) and data has > 800 rows to gain speed from repeated function calls. +- Use a vectorized optimization to improve speed of optimization required for deviance calculation. + # extras 0.7.3.9002 - Remove dependency on MASS package so minimum R version can be brought down to 4.0.0 from 4.3.0. diff --git a/R/dev.R b/R/dev.R index 715a9613..4c2d17b1 100644 --- a/R/dev.R +++ b/R/dev.R @@ -1,3 +1,59 @@ +# Function to optimize a vector of parameter values at the same time +parallel_optimize <- function(f, interval, N, tol = .Machine$double.eps^0.5) { + if (N == 0) { + return(numeric(0)) + } + + lower <- rep(interval[1], N) + upper <- rep(interval[2], N) + threshold <- (upper - lower) * tol + + # Golden ratio + gr <- (sqrt(5) - 1) / 2 + + # Initial points + x1 <- lower + (1 - gr) * (upper - lower) + x2 <- lower + gr * (upper - lower) + f1 <- f(x1) + f2 <- f(x2) + stopifnot(!anyNA(f1), !anyNA(f2)) + + # Iteratively narrow the search interval + while (all((abs(upper - lower) > threshold) | (abs(f2 - f1) > pmax(abs(f1) * tol, 1e-150)))) { + pos_smaller <- f1 < f2 + idx_smaller <- which(pos_smaller) + idx_larger <- which(!pos_smaller) + + upper[idx_smaller] <- x2[idx_smaller] + x2[idx_smaller] <- x1[idx_smaller] + f2[idx_smaller] <- f1[idx_smaller] + x1[idx_smaller] <- lower[idx_smaller] + (1 - gr) * (upper[idx_smaller] - lower[idx_smaller]) + + lower[idx_larger] <- x1[idx_larger] + x1[idx_larger] <- x2[idx_larger] + f1[idx_larger] <- f2[idx_larger] + x2[idx_larger] <- lower[idx_larger] + gr * (upper[idx_larger] - lower[idx_larger]) + + x_new <- ifelse(pos_smaller, x1, x2) + f_new <- f(x_new) + f1[idx_smaller] <- f_new[idx_smaller] + f2[idx_larger] <- f_new[idx_larger] + } + + x_inside <- x1 + x_inside[idx_larger] <- x2[idx_larger] + f_inside <- f1 + f_inside[idx_larger] <- f2[idx_larger] + + x_inside <- ifelse(pos_smaller, x1, x2) + f_inside <- ifelse(pos_smaller, f1, f2) + + x_outside <- ifelse(pos_smaller, lower, upper) + f_outside <- f(x_outside) + + return(ifelse(f_inside < f_outside, x_inside, x_outside)) +} + #' Beta-Binomial Deviances #' #' This parameterization of the beta-binomial distribution uses an expected @@ -21,44 +77,59 @@ #' @examples #' dev_beta_binom(c(0, 1, 2), 10, 0.5, 0.1) dev_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0, res = FALSE) { - opt_beta_binom <- function(prob, x, size = size, theta = theta) { - -log_lik_beta_binom(x = x, size = size, prob = prob, theta = theta) - } + force(x) + args_not_na <- !is.na(x + size + theta) if (length(size) == 1) { - size <- rep(size, length(x)) + size_vec <- size + size_rep <- rep(size, length(x)) + } else { + size_vec <- size[args_not_na] + size_rep <- size } if (length(theta) == 1) { - theta <- rep(theta, length(x)) - } - opt_p <- rep(NA, length(x)) - bol <- !is.na(x) & !is.na(size) & !is.na(theta) - for (i in seq_along(x)) { - if (bol[i] && size[i] < x[i]) { - opt_p[i] <- 1 - } else if (bol[i] && !is.na(bol[i])) { - opt_p[i] <- stats::optimize( - opt_beta_binom, - interval = c(0, 1), - x = x[i], - size = size[i], - theta = theta[i], - tol = 1e-8 - )$minimum - } + theta_vec <- theta + theta_rep <- rep(theta, length(x)) + } else { + theta_vec <- theta[args_not_na] + theta_rep <- theta } + opt_p <- rep(NA, length(args_not_na)) + opt_p[args_not_na] <- parallel_optimize( + f = make_opt_beta_binom(x[args_not_na], size_vec, theta_vec), + interval = c(0, 1), + N = sum(args_not_na) + ) dev1 <- log_lik_beta_binom(x = x, size = size, prob = opt_p, theta = theta) dev2 <- log_lik_beta_binom(x = x, size = size, prob = prob, theta = theta) dev <- dev1 - dev2 dev[dev < 0 & dev > -1e-7] <- 0 dev <- dev * 2 - use_binom <- (!is.na(theta) & theta == 0) | - (!is.na(x) & !is.na(size) & x == 0 & size == 0) + use_binom <- (!is.na(theta_rep) & theta_rep == 0) | + (!is.na(x) & !is.na(size_rep) & x == 0 & size_rep == 0) dev_binom <- dev_binom(x = x, size = size, prob = prob, res = FALSE) dev[use_binom] <- dev_binom[use_binom] if (vld_false(res)) { return(dev) } - dev_res(x, ifelse(use_binom, size * prob, size * opt_p), dev) + dev_res(x, size_rep * prob, dev) +} + +# Objective function for optimization. +make_opt_beta_binom <- function(x, size, theta) { + force(x) + force(size) + force(theta) + function(prob) { + out <- -log_lik_beta_binom( + x = x, + prob = prob, + size = size, + theta = theta, + memoize = TRUE + ) + out[is.nan(out)] <- Inf + out + } } #' Bernoulli Deviances diff --git a/R/log-lik.R b/R/log-lik.R index 4d84dbb2..60d1af27 100644 --- a/R/log-lik.R +++ b/R/log-lik.R @@ -12,6 +12,7 @@ #' #' @inheritParams params #' @param x A non-negative whole numeric vector of values. +#' @param memoize Whether or not to memoize the function. #' #' @return An numeric vector of the corresponding log-likelihoods. #' @family log_lik_dist @@ -19,28 +20,75 @@ #' #' @examples #' log_lik_beta_binom(c(0, 1, 2), 3, 0.5, 0) -log_lik_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0) { +log_lik_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0, memoize = FALSE) { alpha <- prob * 2 * (1 / theta) beta <- (1 - prob) * 2 * (1 / theta) - lbeta_binom <- lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) + + + # Memoise use case: + # Posterior_predictive_check calls this function repeatedly with + # x and size unchanged; memoize it to reduce repeated calls + # when length(x) is large enough to outweigh the overhead required by memoize. + # For length(x) < 800, memoize is slower + # See detail here https://github.com/poissonconsulting/extras/issues/63 + if (memoize && length(x) >= 800) { + lgamma_size_x <- lgamma_size_x(size, x) + } else { + lgamma_size_x <- lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) + } + lbeta_binom <- lgamma_size_x + lgamma(x + alpha) + lgamma(size - x + beta) - lgamma(size + alpha + beta) + lgamma(alpha + beta) - lgamma(alpha) - lgamma(beta) - bol <- !is.na(x + size + prob + theta) - lbeta_binom[bol & ((x == 0 & prob == 0) | (x == size & prob == 1))] <- 0 - lbeta_binom[bol & x != 0 & prob == 0] <- -Inf - lbeta_binom[bol & x != size & prob == 1] <- -Inf - lbeta_binom[bol & x > size] <- -Inf - bol_theta <- !is.na(theta) - lbeta_binom[bol_theta & theta < 0] <- NaN - use_binom <- bol_theta & theta == 0 - if (any(use_binom)) { - lbinom <- log_lik_binom(x, size = size, prob = prob) - lbeta_binom[use_binom] <- lbinom[use_binom] - } - if (length(bol) == 0) { - lbeta_binom <- numeric(0) + + args_na <- is.na(x + size + prob + theta) + length_args_na <- length(args_na) + if (length_args_na == 1) { + if (args_na) { + return(NA_real_) + } + if (prob == 0) { + if (x == 0) { + lbeta_binom <- 0 + } else { + lbeta_binom <- -Inf + } + } else if (prob == 1) { + if (x == size) { + lbeta_binom <- 0 + } else { + lbeta_binom <- -Inf + } + } + if (x > size) { + lbeta_binom <- -Inf + } + if (theta == 0) { + lbeta_binom <- log_lik_binom(x = x, size = size, prob = prob) + } else if (theta < 0) { + lbeta_binom <- NaN + } + lbeta_binom + } else if (length_args_na == 0) { + numeric(0) + } else { + args_not_na <- !args_na + lbeta_binom[args_not_na & ((x == 0 & prob == 0) | (x == size & prob == 1))] <- 0 + lbeta_binom[args_not_na & x != 0 & prob == 0] <- -Inf + lbeta_binom[args_not_na & x != size & prob == 1] <- -Inf + lbeta_binom[args_not_na & x > size] <- -Inf + theta_not_na <- !is.na(theta) + lbeta_binom[theta_not_na & theta < 0] <- NaN + use_binom <- theta_not_na & theta == 0 + if (any(use_binom)) { + lbinom <- log_lik_binom(x, size = size, prob = prob) + lbeta_binom[use_binom] <- lbinom[use_binom] + } + lbeta_binom } - lbeta_binom +} + +# Function to memoize (called repeatedly for non-changing values of size and x) +lgamma_size_x <- function(size, x) { + lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) } #' Bernoulli Log-Likelihood diff --git a/R/zzz.R b/R/zzz.R new file mode 100644 index 00000000..9e4d08ad --- /dev/null +++ b/R/zzz.R @@ -0,0 +1,5 @@ +.onLoad <- function(libname, pkgname) { + if (requireNamespace("memoise", quietly = TRUE)) { + lgamma_size_x <<- memoise::memoise(lgamma_size_x) + } +} diff --git a/README.Rmd b/README.Rmd index 5bca6aa9..2dd2cd7c 100644 --- a/README.Rmd +++ b/README.Rmd @@ -13,6 +13,7 @@ knitr::opts_chunk$set( ) ``` + # extras extras website @@ -32,19 +33,6 @@ calculate deviance residuals as well as R translations of some BUGS (Bayesian Using Gibbs Sampling), JAGS (Just Another Gibbs Sampler), STAN and TMB (Template Model Builder) functions. -## Installation - - -```{r, eval=FALSE, echo=FALSE} -install.packages("extras") -``` - -To install the developmental version from [GitHub](https://github.com/poissonconsulting/extras) -```{r, eval=FALSE} -# install.packages("pak") -pak::pak("poissonconsulting/extras") -``` - ## Demonstration ### Summarise MCMC Samples @@ -67,6 +55,7 @@ svalue(rnorm(100, mean = 3)) Implemented distributions with functions to draw random samples, calculate log-likelihoods, and calculate deviance residuals for include: - Bernoulli +- Binomial - Beta-binomial - Gamma - Gamma-Poisson @@ -106,6 +95,37 @@ numericise( ) ``` +## Installation + + +## Information + +For more information see the [Get Started](https://poissonconsulting.github.io/chk/articles/chk.html) vignette. + +## Installation + +### Release + +To install the release version from [CRAN](https://CRAN.R-project.org/package=extras). +```r +install.packages("extras") +``` + +The website for the release version is at . + +### Development + +To install the development version from [GitHub](https://github.com/poissonconsulting/extras) +```r +# install.packages("remotes") +remotes::install_github("poissonconsulting/extras") +``` + +or from [r-universe](https://poissonconsulting.r-universe.dev/extras). +```r +install.packages("extras", repos = c("https://poissonconsulting.r-universe.dev", "https://cloud.r-project.org")) +``` + ## References Greenland, S. 2019. Valid P-Values Behave Exactly as They Should: Some Misleading Criticisms of P-Values and Their Resolution With S-Values. The American Statistician 73(sup1): 106–114. @@ -120,4 +140,3 @@ Please report any [issues](https://github.com/poissonconsulting/extras/issues). Please note that the extras project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms. - diff --git a/README.md b/README.md index e95b0084..a3325974 100644 --- a/README.md +++ b/README.md @@ -23,18 +23,6 @@ distributions and calculate deviance residuals as well as R translations of some BUGS (Bayesian Using Gibbs Sampling), JAGS (Just Another Gibbs Sampler), STAN and TMB (Template Model Builder) functions. -## Installation - - - -To install the developmental version from -[GitHub](https://github.com/poissonconsulting/extras) - -``` r -# install.packages("pak") -pak::pak("poissonconsulting/extras") -``` - ## Demonstration ### Summarise MCMC Samples @@ -68,6 +56,7 @@ Implemented distributions with functions to draw random samples, calculate log-likelihoods, and calculate deviance residuals for include: - Bernoulli +- Binomial - Beta-binomial - Gamma - Gamma-Poisson @@ -115,6 +104,44 @@ numericise( #> [2,] 0 2 10958 61 ``` +## Installation + +## Information + +For more information see the [Get +Started](https://poissonconsulting.github.io/chk/articles/chk.html) +vignette. + +## Installation + +### Release + +To install the release version from +[CRAN](https://CRAN.R-project.org/package=extras). + +``` r +install.packages("extras") +``` + +The website for the release version is at +. + +### Development + +To install the development version from +[GitHub](https://github.com/poissonconsulting/extras) + +``` r +# install.packages("remotes") +remotes::install_github("poissonconsulting/extras") +``` + +or from [r-universe](https://poissonconsulting.r-universe.dev/extras). + +``` r +install.packages("extras", repos = c("https://poissonconsulting.r-universe.dev", "https://cloud.r-project.org")) +``` + ## References Greenland, S. 2019. Valid P-Values Behave Exactly as They Should: Some diff --git a/_pkgdown.yml b/_pkgdown.yml index 70859e68..572db40d 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,5 +1,11 @@ -destination: docs url: https://poissonconsulting.github.io/extras/ + +template: + package: poissontemplate + +development: + mode: auto + reference: - title: Manipulate Objects desc: Functions to numericise and fill objects @@ -187,5 +193,4 @@ reference: - '`pextreme`' - '`sextreme`' - '`as_list_unnamed`' -template: - package: poissontemplate + diff --git a/extras.Rproj b/extras.Rproj index 88ff2b5d..9e36e1db 100644 --- a/extras.Rproj +++ b/extras.Rproj @@ -1,4 +1,5 @@ Version: 1.0 +ProjectId: 8629deed-40ee-4c36-9796-bc4f821454de RestoreWorkspace: Default SaveWorkspace: Default diff --git a/man/dev_skewnorm.Rd b/man/dev_skewnorm.Rd index 9983b41c..1726fe2d 100644 --- a/man/dev_skewnorm.Rd +++ b/man/dev_skewnorm.Rd @@ -25,7 +25,7 @@ An numeric vector of the corresponding deviances or deviance residuals. Skew Normal Deviances } \examples{ -\dontshow{if (rlang::is_installed("sn")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("sn")) withAutoprint(\{ # examplesIf} dev_skewnorm(c(-2:2)) dev_skewnorm(-2:2, 0, 1, 5) dev_skewnorm(-2:2, 0, 1, 5, res = TRUE) diff --git a/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd index e55cbea6..880de3e4 100644 --- a/man/log_lik_beta_binom.Rd +++ b/man/log_lik_beta_binom.Rd @@ -4,7 +4,7 @@ \alias{log_lik_beta_binom} \title{Beta-Binomial Log-Likelihood} \usage{ -log_lik_beta_binom(x, size = 1, prob = 0.5, theta = 0) +log_lik_beta_binom(x, size = 1, prob = 0.5, theta = 0, memoize = FALSE) } \arguments{ \item{x}{A non-negative whole numeric vector of values.} @@ -16,6 +16,8 @@ success.} \item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).} + +\item{memoize}{Whether or not to memoize the function.} } \value{ An numeric vector of the corresponding log-likelihoods. diff --git a/man/log_lik_skewnorm.Rd b/man/log_lik_skewnorm.Rd index ea6eae50..f468eb28 100644 --- a/man/log_lik_skewnorm.Rd +++ b/man/log_lik_skewnorm.Rd @@ -22,7 +22,7 @@ An numeric vector of the corresponding log-likelihoods. Skew Normal Log-Likelihood } \examples{ -\dontshow{if (rlang::is_installed("sn")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("sn")) withAutoprint(\{ # examplesIf} log_lik_skewnorm(c(-2:2)) log_lik_skewnorm(c(-2:2), shape = -2) log_lik_skewnorm(c(-2:2), shape = 2) diff --git a/man/numericise.Rd b/man/numericise.Rd index 1d7215b3..f500d844 100644 --- a/man/numericise.Rd +++ b/man/numericise.Rd @@ -100,7 +100,7 @@ numericise(as.Date("1972-01-01")) numericise(as.POSIXct("1972-01-01", tz = "UTC")) # hms -\dontshow{if (rlang::is_installed("hms")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("hms")) withAutoprint(\{ # examplesIf} numericise(hms::as_hms("00:01:03")) \dontshow{\}) # examplesIf} diff --git a/man/ran_skewnorm.Rd b/man/ran_skewnorm.Rd index b20f753e..805c4d11 100644 --- a/man/ran_skewnorm.Rd +++ b/man/ran_skewnorm.Rd @@ -23,7 +23,7 @@ A numeric vector of the random samples. Skew Normal Random Samples } \examples{ -\dontshow{if (rlang::is_installed("sn")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("sn")) withAutoprint(\{ # examplesIf} ran_skewnorm(10, shape = -1) ran_skewnorm(10, shape = 0) ran_skewnorm(10, shape = 1) diff --git a/man/res_skewnorm.Rd b/man/res_skewnorm.Rd index ecee901e..0d2c53c5 100644 --- a/man/res_skewnorm.Rd +++ b/man/res_skewnorm.Rd @@ -27,7 +27,7 @@ An numeric vector of the corresponding residuals. Skew Normal Residuals } \examples{ -\dontshow{if (rlang::is_installed("sn")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("sn")) withAutoprint(\{ # examplesIf} res_skewnorm(c(-2:2)) \dontshow{\}) # examplesIf} } diff --git a/man/sens_skewnorm.Rd b/man/sens_skewnorm.Rd index 262637a7..fb5605d2 100644 --- a/man/sens_skewnorm.Rd +++ b/man/sens_skewnorm.Rd @@ -24,7 +24,7 @@ Expands (\code{sd_mult > 1}) or reduces (\code{sd_mult < 1}) the standard deviat of the Skew Normal distribution without changing the mean. } \examples{ -\dontshow{if (rlang::is_installed("sn")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("sn")) withAutoprint(\{ # examplesIf} sens_skewnorm(10, 3, -1, 2) sens_skewnorm(10, 3, 3, 0.8) \dontshow{\}) # examplesIf} diff --git a/man/skewnorm.Rd b/man/skewnorm.Rd index 77747585..59007b5c 100644 --- a/man/skewnorm.Rd +++ b/man/skewnorm.Rd @@ -41,7 +41,7 @@ generate.} Skew-Normal Distribution } \examples{ -\dontshow{if (rlang::is_installed("sn")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (rlang::is_installed("sn")) withAutoprint(\{ # examplesIf} dskewnorm(x = -2:2, mean = 0, sd = 1, shape = 0.1) dskewnorm(x = -2:2, mean = 0, sd = 1, shape = -1) qskewnorm(p = c(0.1, 0.4), mean = 0, sd = 1, shape = 0.1) diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R index 19d52086..6bd2a5e0 100644 --- a/tests/testthat/test-dev.R +++ b/tests/testthat/test-dev.R @@ -12,8 +12,8 @@ test_that("beta_binom known values", { expect_equal(dev_beta_binom(1, 1, 0, 0.5), Inf) expect_equal(dev_beta_binom(1, 1, 0.5, 0), 1.38629436111989) expect_identical(dev_beta_binom(1, 1, 0.5, 0), dev_binom(1, 1, 0.5)) - expect_equal(dev_beta_binom(1, 1, 0.5, 1), 1.38629430121326) - expect_equal(dev_beta_binom(1, 1, 0.5, 0.5), 1.38629430121326) + expect_equal(dev_beta_binom(1, 1, 0.5, 1), 1.38629436111989) + expect_equal(dev_beta_binom(1, 1, 0.5, 0.5), 1.38629436111989) expect_equal(dev_beta_binom(1, 2, 0.2, 0), 0.892574205256839) expect_equal(dev_beta_binom(1, 2, 0.2, 1), 0.892574205256839) expect_equal(dev_beta_binom(1, 2, 0.2, 0.5), 0.892574205256839) @@ -39,8 +39,8 @@ test_that("beta_binom known values", { expect_equal(dev_beta_binom(0, 2, 0.5, 0.1), 2.67954869096997) expect_equal(dev_beta_binom(0, 2, 0.5, 0.5), 2.40794560865187) expect_equal(dev_beta_binom(0, 2, 0.1), 0.421442062631305) - expect_equal(dev_beta_binom(0, 2, 0.1, 0.1), 0.410887933834857) - expect_equal(dev_beta_binom(0, 2, 0.1, 0.5), 0.377484235738089) + expect_equal(dev_beta_binom(0, 2, 0.1, 0.1), 0.410887948429604) + expect_equal(dev_beta_binom(0, 2, 0.1, 0.5), 0.377484249193745) expect_equal(dev_beta_binom(1, 2, 1, 10), Inf) expect_equal(dev_beta_binom(2, 2, 1, 10), 0) expect_equal(dev_beta_binom(0, 2, 0.5, 10), 1.56031711509915) @@ -49,7 +49,7 @@ test_that("beta_binom known values", { test_that("beta_binom vectorized", { expect_equal(dev_beta_binom(0:3, 5, 0, 0), dev_binom(0:3, 5, 0)) - expect_equal(dev_beta_binom(c(0, 1, 3, 0), 3, 0.5, 0.5), c(3.21887580642895, 0.179750127270295, 3.2188756770985, 3.21887580642895)) + expect_equal(dev_beta_binom(c(0, 1, 3, 0), 3, 0.5, 0.5), c(3.2188758248682, 0.179750127270293, 3.2188758248682, 3.2188758248682)) expect_equal(dev_beta_binom(0:3, 0:3, rep(1, 4), 0), rep(0, 4)) expect_equal(dev_beta_binom(0:3, 1:4, seq(0, 1, length.out = 4), 0:3), c(0, 0.235566071312767, 0.076961041136129, Inf)) }) @@ -256,17 +256,16 @@ test_that("gamma_pois log_lik", { # }) test_that("gamma_pois deviance snapshot", { - expect_snapshot( - { - withr::with_seed( - 101, { - x <- ran_gamma_pois(10000, 3, 0.5) - } - ) - deviance <- dev_gamma_pois(x, 3, 0.5) - deviance - } - ) + expect_snapshot({ + withr::with_seed( + 101, + { + x <- ran_gamma_pois(10000, 3, 0.5) + } + ) + deviance <- dev_gamma_pois(x, 3, 0.5) + deviance + }) }) test_that("gamma_pois ran", { diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R index 5419a119..70de5b18 100644 --- a/tests/testthat/test-log-lik.R +++ b/tests/testthat/test-log-lik.R @@ -165,10 +165,57 @@ test_that("beta_binom known values", { expect_equal(log_lik_beta_binom(1, 2, 0.2, 1), -1.54489939129653) expect_equal(log_lik_beta_binom(2, 2, 0.2, 10), -1.75253875607477) expect_equal(log_lik_beta_binom(1, 2, 0.5), -0.693147180559945) + expect_equal(log_lik_beta_binom(10, 2, 0.5, 0.1), -Inf) + expect_equal(log_lik_beta_binom(10, 2, 0.5, -0.1), NaN) expect_equal(log_lik_beta_binom(1, 2, 0.5), log_lik_binom(1, 2, 0.5)) expect_equal(log_lik_beta_binom(1, 2, 0.2), dbinom(1, 2, 0.2, log = TRUE)) }) +test_that("beta binomial deviance function is memoized", { + expect_true(memoise::is.memoized(lgamma_size_x)) +}) + +test_that("lgamma_size_x produces expected outputs", { + size <- 100 + x <- 1 + expect_equal( + lgamma_size_x(size, x), + lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) + ) + size <- 100 + x <- 1:100 + expect_equal( + lgamma_size_x(size, x), + lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) + ) + size <- 201:300 + x <- 1:100 + expect_equal( + lgamma_size_x(size, x), + lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) + ) +}) + +test_that("beta_binom memoized function gives same outputs as non-memoized function", { + expect_equal( + log_lik_beta_binom(1:100, 200, 0.5, 0.1, memoize = FALSE), + log_lik_beta_binom(1:100, 200, 0.5, 0.1, memoize = TRUE) + ) + expect_equal( + log_lik_beta_binom(1:100, 200, 0.1, 0, memoize = FALSE), + log_lik_beta_binom(1:100, 200, 0.1, 0, memoize = TRUE) + ) + withr::with_seed( + 101, + { + x <- ran_beta_binom(1e7, 10, 0.5, 0.1) + } + ) + time_no_mem <- system.time(log_lik_beta_binom(x, 10, 0.5, 0.2, memoize = FALSE))[3] + time_mem <- system.time(log_lik_beta_binom(x, 10, 0.5, 0.2, memoize = TRUE))[3] + expect_true(time_mem < time_no_mem) +}) + test_that("beta_binom vectorized", { expect_equal(log_lik_beta_binom(0:3, 2, 0.5, 0), c(-1.38629436111989, -0.693147180559945, -1.38629436111989, -Inf)) expect_equal(log_lik_beta_binom(0:3, 2, 0.5, 0), log_lik_binom(0:3, 2, 0.5)) diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 1bb8da5f..0bbf882c 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -509,9 +509,9 @@ test_that("res_beta_binom", { expect_equal(res_beta_binom(0, 1, 0.7), -1.55175565365552) expect_equal(res_beta_binom(1, 1, 0.7), 0.844600430900592) expect_equal(res_beta_binom(0, 2, 0.5, 0.5), -1.55175564931989) - expect_equal(res_beta_binom(1, 2, 0.5, 0.5), 2.1073424255447e-08) + expect_equal(res_beta_binom(1, 2, 0.5, 0.5), 0) expect_equal(res_beta_binom(0, 2, 0.7, 0.5), -2.01243800222306) - expect_equal(res_beta_binom(1, 2, 0.7, 0.5), 0.590513991612016) + expect_equal(res_beta_binom(1, 2, 0.7, 0.5), -0.590513991612016) expect_identical(res_beta_binom(1, 2, 0.5), 0) expect_identical(res_beta_binom(5, 10, 0.5), 0) expect_identical(res_beta_binom(1, 2, 0.5, 2), 0) @@ -533,7 +533,7 @@ test_that("res_beta_binom", { ) ) expect_equal(res_beta_binom(0, 2, 0.5, type = "dev"), -1.6651092223154) - expect_equal(res_beta_binom(0, 2, 0.5, 10, type = "dev"), -1.24912653737637) + expect_equal(res_beta_binom(0, 2, 0.5, 10, type = "dev"), -1.24912654086732) expect_equal(res_beta_binom(0, 2, 0.5, type = "raw"), -1) expect_equal(res_beta_binom(0, 2, 0.5, 10, type = "raw"), -1) set.seed(101)