From a386e0e6872f6e4343bee620c064fc4e673d455e Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Thu, 8 Aug 2024 11:21:46 -0700
Subject: [PATCH 01/20] Add binomial distribution to README

---
 README.Rmd | 1 +
 README.md  | 1 +
 2 files changed, 2 insertions(+)

diff --git a/README.Rmd b/README.Rmd
index b9b264ca..e82eb373 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -67,6 +67,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
diff --git a/README.md b/README.md
index 07dcca36..403bea85 100644
--- a/README.md
+++ b/README.md
@@ -68,6 +68,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

From 3ca7c00a5351638af5189b147e8027b819f20db5 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Thu, 8 Aug 2024 11:23:16 -0700
Subject: [PATCH 02/20] Optimize log-likelihood function for beta-binomial
 deviance calculation #63

---
 R/log-lik.R | 60 +++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 45 insertions(+), 15 deletions(-)

diff --git a/R/log-lik.R b/R/log-lik.R
index e7c138ea..825c42e8 100644
--- a/R/log-lik.R
+++ b/R/log-lik.R
@@ -25,22 +25,52 @@ log_lik_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0) {
   lbeta_binom <- lgamma(size + 1) - lgamma(x + 1) - lgamma(size - x + 1) +
     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
 }
 
 #' Bernoulli Log-Likelihood

From 8632d1dae0cf71d773d2735635a8493e5e98b257 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Thu, 8 Aug 2024 11:24:33 -0700
Subject: [PATCH 03/20] Styler

---
 tests/testthat/test-chk-index.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/testthat/test-chk-index.R b/tests/testthat/test-chk-index.R
index 10475ae6..e441704c 100644
--- a/tests/testthat/test-chk-index.R
+++ b/tests/testthat/test-chk-index.R
@@ -15,14 +15,14 @@ test_that("chk_index", {
   )
 })
 
-test_that("chk_index errors with x = NA",{
+test_that("chk_index errors with x = NA", {
   expect_snapshot(
     error = TRUE,
     chk_index(NA)
   )
 })
 
-test_that("chk_index errors with empty x",{
+test_that("chk_index errors with empty x", {
   expect_snapshot(
     error = TRUE,
     chk_index(integer(0))

From 0a51424579bfdc26fc5d49a27f2536175c81e651 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 12 Aug 2024 09:29:28 -0700
Subject: [PATCH 04/20] Add memoized function to repeated portion of
 log-likelihood

---
 DESCRIPTION               |  1 +
 R/log-lik.R               | 22 ++++++++++++++++++++--
 R/zzz.R                   |  5 +++++
 man/log_lik_beta_binom.Rd |  4 +++-
 4 files changed, 29 insertions(+), 3 deletions(-)
 create mode 100644 R/zzz.R

diff --git a/DESCRIPTION b/DESCRIPTION
index ee5e566c..928f5b4c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -34,6 +34,7 @@ Suggests:
     ggplot2,
     knitr,
     MASS,
+    memoise,
     rlang,
     rmarkdown,
     scales,
diff --git a/R/log-lik.R b/R/log-lik.R
index 825c42e8..c1c782c5 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,10 +20,22 @@
 #'
 #' @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
+  # https://poissonconsulting.slack.com/archives/C07FC57Q346/p1723248220478269
+  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)
 
@@ -73,6 +86,11 @@ log_lik_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0) {
   }
 }
 
+# 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
 #'
 #' @inheritParams params
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/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd
index e84a44f9..761d9e84 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.}
@@ -14,6 +14,8 @@ log_lik_beta_binom(x, size = 1, prob = 0.5, theta = 0)
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of 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.

From df94a1b77553ca9446a4be560043e31c6a5d9338 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 12 Aug 2024 09:30:19 -0700
Subject: [PATCH 05/20] - Use vectorized optimization function to speed up
 deviance calculation for the beta-binomial distribution

---
 R/dev.R | 160 +++++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 136 insertions(+), 24 deletions(-)

diff --git a/R/dev.R b/R/dev.R
index 0bcf9457..9a6389d1 100644
--- a/R/dev.R
+++ b/R/dev.R
@@ -1,3 +1,100 @@
+# Define a custom optimize function in R
+optimize_R <-
+  custom_optimize <- function(f, interval, tol = .Machine$double.eps^0.25, maximum = FALSE) {
+    lower <- interval[1]
+    upper <- interval[2]
+
+    # 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)
+
+    # Iteratively narrow the search interval
+    while (abs(upper - lower) > tol * (abs(x1) + abs(x2))) {
+      if ((f1 < f2) != maximum) {
+        upper <- x2
+        x2 <- x1
+        f2 <- f1
+        x1 <- lower + (1 - gr) * (upper - lower)
+        f1 <- f(x1)
+      } else {
+        lower <- x1
+        x1 <- x2
+        f1 <- f2
+        x2 <- lower + gr * (upper - lower)
+        f2 <- f(x2)
+      }
+    }
+
+    # Return the minimum (or maximum) point and the function value at that point
+    if (f1 < f2) {
+      return(list(minimum = x1, objective = f1))
+    } else {
+      return(list(minimum = x2, objective = f2))
+    }
+  }
+
+# Function to optimize a vector of parameter values at the same time
+# Modified from `optimize_R()` above.
+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 +118,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

From dfd14061f85dd61ee04ff25d3d19bf4f83edf4b5 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 12 Aug 2024 09:32:41 -0700
Subject: [PATCH 06/20] Improved optimization has slightly different values -
 update tests accordingly

---
 tests/testthat/test-dev.R | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R
index 062c4d71..44269c5d 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))
 })

From ecfeebc57d958c761e2ba41ee16279d4117a1d5e Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 12 Aug 2024 09:32:53 -0700
Subject: [PATCH 07/20] Improved optimization has slightly different values -
 update tests accordingly

---
 tests/testthat/test-res.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R
index 44552452..9a2672d6 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)

From 3455738524823042ab1aed31928ca1cfa42a235d Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Tue, 13 Aug 2024 09:23:37 -0700
Subject: [PATCH 08/20] Comment out original R optimize function (not used)

---
 R/dev.R | 78 ++++++++++++++++++++++++++++-----------------------------
 1 file changed, 39 insertions(+), 39 deletions(-)

diff --git a/R/dev.R b/R/dev.R
index 9a6389d1..dc612c13 100644
--- a/R/dev.R
+++ b/R/dev.R
@@ -1,42 +1,42 @@
-# Define a custom optimize function in R
-optimize_R <-
-  custom_optimize <- function(f, interval, tol = .Machine$double.eps^0.25, maximum = FALSE) {
-    lower <- interval[1]
-    upper <- interval[2]
-
-    # 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)
-
-    # Iteratively narrow the search interval
-    while (abs(upper - lower) > tol * (abs(x1) + abs(x2))) {
-      if ((f1 < f2) != maximum) {
-        upper <- x2
-        x2 <- x1
-        f2 <- f1
-        x1 <- lower + (1 - gr) * (upper - lower)
-        f1 <- f(x1)
-      } else {
-        lower <- x1
-        x1 <- x2
-        f1 <- f2
-        x2 <- lower + gr * (upper - lower)
-        f2 <- f(x2)
-      }
-    }
-
-    # Return the minimum (or maximum) point and the function value at that point
-    if (f1 < f2) {
-      return(list(minimum = x1, objective = f1))
-    } else {
-      return(list(minimum = x2, objective = f2))
-    }
-  }
+# # Define a custom optimize function in R
+# optimize_R <-
+#   custom_optimize <- function(f, interval, tol = .Machine$double.eps^0.25, maximum = FALSE) {
+#     lower <- interval[1]
+#     upper <- interval[2]
+#
+#     # 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)
+#
+#     # Iteratively narrow the search interval
+#     while (abs(upper - lower) > tol * (abs(x1) + abs(x2))) {
+#       if ((f1 < f2) != maximum) {
+#         upper <- x2
+#         x2 <- x1
+#         f2 <- f1
+#         x1 <- lower + (1 - gr) * (upper - lower)
+#         f1 <- f(x1)
+#       } else {
+#         lower <- x1
+#         x1 <- x2
+#         f1 <- f2
+#         x2 <- lower + gr * (upper - lower)
+#         f2 <- f(x2)
+#       }
+#     }
+#
+#     # Return the minimum (or maximum) point and the function value at that point
+#     if (f1 < f2) {
+#       return(list(minimum = x1, objective = f1))
+#     } else {
+#       return(list(minimum = x2, objective = f2))
+#     }
+#   }
 
 # Function to optimize a vector of parameter values at the same time
 # Modified from `optimize_R()` above.

From 6a0284b59b5579e5073fb34affa29716917a95cd Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Tue, 13 Aug 2024 09:24:45 -0700
Subject: [PATCH 09/20] Add tests for memoized function and two cases missing
 test coverage in scalar log likelihood calculation

---
 tests/testthat/test-log-lik.R | 46 +++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)

diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R
index 824bc6f6..bd7b2595 100644
--- a/tests/testthat/test-log-lik.R
+++ b/tests/testthat/test-log-lik.R
@@ -165,10 +165,56 @@ 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))

From e789cc49e8c8af612a7c66d5daa73ad3198b4dd9 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Tue, 13 Aug 2024 14:39:18 -0700
Subject: [PATCH 10/20] Move `optimize_R()` function to the scripts folder for
 reference

---
 R/dev.R | 41 -----------------------------------------
 1 file changed, 41 deletions(-)

diff --git a/R/dev.R b/R/dev.R
index dc612c13..22580890 100644
--- a/R/dev.R
+++ b/R/dev.R
@@ -1,45 +1,4 @@
-# # Define a custom optimize function in R
-# optimize_R <-
-#   custom_optimize <- function(f, interval, tol = .Machine$double.eps^0.25, maximum = FALSE) {
-#     lower <- interval[1]
-#     upper <- interval[2]
-#
-#     # 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)
-#
-#     # Iteratively narrow the search interval
-#     while (abs(upper - lower) > tol * (abs(x1) + abs(x2))) {
-#       if ((f1 < f2) != maximum) {
-#         upper <- x2
-#         x2 <- x1
-#         f2 <- f1
-#         x1 <- lower + (1 - gr) * (upper - lower)
-#         f1 <- f(x1)
-#       } else {
-#         lower <- x1
-#         x1 <- x2
-#         f1 <- f2
-#         x2 <- lower + gr * (upper - lower)
-#         f2 <- f(x2)
-#       }
-#     }
-#
-#     # Return the minimum (or maximum) point and the function value at that point
-#     if (f1 < f2) {
-#       return(list(minimum = x1, objective = f1))
-#     } else {
-#       return(list(minimum = x2, objective = f2))
-#     }
-#   }
-
 # Function to optimize a vector of parameter values at the same time
-# Modified from `optimize_R()` above.
 parallel_optimize <- function(f, interval, N, tol = .Machine$double.eps^0.5) {
   if (N == 0) {
     return(numeric(0))

From 17386a4bb0dd267516e316fbee430d59bac37c27 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Tue, 13 Aug 2024 14:42:27 -0700
Subject: [PATCH 11/20] Remove link

---
 R/log-lik.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/log-lik.R b/R/log-lik.R
index c1c782c5..15692b10 100644
--- a/R/log-lik.R
+++ b/R/log-lik.R
@@ -29,7 +29,7 @@ log_lik_beta_binom <- function(x, size = 1, prob = 0.5, theta = 0, memoize = FAL
   # 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
-  # https://poissonconsulting.slack.com/archives/C07FC57Q346/p1723248220478269
+  # See detail here https://github.com/poissonconsulting/extras/issues/63
   if (memoize && length(x) >= 800) {
     lgamma_size_x <- lgamma_size_x(size, x)
   } else {

From 4aef63eae63a93efec0e34e2955e8776a4784961 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Tue, 17 Sep 2024 13:10:43 -0700
Subject: [PATCH 12/20] Update documentation

---
 man/log_lik_beta_binom.Rd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd
index 31573ab0..880de3e4 100644
--- a/man/log_lik_beta_binom.Rd
+++ b/man/log_lik_beta_binom.Rd
@@ -16,7 +16,7 @@ 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{

From 1bd7cd811fd99e04ba000451c2c2577199a1aec5 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 15:08:13 -0800
Subject: [PATCH 13/20] add ProjectId

---
 extras.Rproj | 1 +
 1 file changed, 1 insertion(+)

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

From f27f8abd1568d315938a89d7bb071120e90a1429 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 15:22:11 -0800
Subject: [PATCH 14/20] Roxygen to 7.3.2.9000 and update associated
 documentation

---
 DESCRIPTION             | 2 +-
 man/dev_skewnorm.Rd     | 2 +-
 man/log_lik_skewnorm.Rd | 2 +-
 man/numericise.Rd       | 2 +-
 man/ran_skewnorm.Rd     | 2 +-
 man/res_skewnorm.Rd     | 2 +-
 man/sens_skewnorm.Rd    | 2 +-
 man/skewnorm.Rd         | 2 +-
 8 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 1e8a0e22..db67728a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -46,6 +46,6 @@ Config/testthat/edition: 3
 Encoding: UTF-8
 Language: en-US
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.2.9000
 VignetteBuilder: knitr
 Config/Needs/website: poissonconsulting/poissontemplate
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_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)

From 90ffbf70d9db56cf2716453c71d937b4bec54015 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 15:31:13 -0800
Subject: [PATCH 15/20] fledge: Bump version to 0.8.0

---
 DESCRIPTION |  2 +-
 NEWS.md     | 16 ++++++++++++++++
 2 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index db67728a..53b611b2 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")),
diff --git a/NEWS.md b/NEWS.md
index 5fb4d984..61f3ccab 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,21 @@
 <!-- NEWS.md is maintained by https://fledge.cynkra.com, contributors should not edit this file -->
 
+# extras 0.8.0
+
+- Up-to-date with main branch.
+
+  Merge remote-tracking branch 'origin/main' into f-optimize-betabin
+
+  # Conflicts:
+  # DESCRIPTION
+
+- Merge branch 'main' into f-optimize-betabin.
+
+- Use vectorized optimization function to speed up deviance calculation for the beta-binomial distribution.
+
+- Optimize log-likelihood function for beta-binomial 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.

From 5b60d948511e8579412d022a4abe60aa2371f4c8 Mon Sep 17 00:00:00 2001
From: Nicole Hill <nicole@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 15:34:27 -0800
Subject: [PATCH 16/20] Update NEWS file

---
 NEWS.md | 16 +++-------------
 1 file changed, 3 insertions(+), 13 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 61f3ccab..992fd516 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,19 +2,9 @@
 
 # extras 0.8.0
 
-- Up-to-date with main branch.
-
-  Merge remote-tracking branch 'origin/main' into f-optimize-betabin
-
-  # Conflicts:
-  # DESCRIPTION
-
-- Merge branch 'main' into f-optimize-betabin.
-
-- Use vectorized optimization function to speed up deviance calculation for the beta-binomial distribution.
-
-- Optimize log-likelihood function for beta-binomial deviance calculation.
-
+- 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
 

From 931c4cc6c53667f9d4b56d35322f1a2ca99d6899 Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 19:14:46 -0800
Subject: [PATCH 17/20] style

---
 tests/testthat/test-dev.R     | 21 ++++++++++-----------
 tests/testthat/test-log-lik.R |  3 ++-
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R
index fdce6d48..6bd2a5e0 100644
--- a/tests/testthat/test-dev.R
+++ b/tests/testthat/test-dev.R
@@ -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 50a67575..70de5b18 100644
--- a/tests/testthat/test-log-lik.R
+++ b/tests/testthat/test-log-lik.R
@@ -206,7 +206,8 @@ test_that("beta_binom memoized function gives same outputs as non-memoized funct
     log_lik_beta_binom(1:100, 200, 0.1, 0, memoize = TRUE)
   )
   withr::with_seed(
-    101, {
+    101,
+    {
       x <- ran_beta_binom(1e7, 10, 0.5, 0.1)
     }
   )

From 80c177fadc8d4cac07ee687808873830eef3ad58 Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 19:31:28 -0800
Subject: [PATCH 18/20] tidy description and rearrange pkgdown

---
 DESCRIPTION  | 5 +++--
 _pkgdown.yml | 8 +++++---
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 53b611b2..332cc5f7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -42,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.9000
-VignetteBuilder: knitr
-Config/Needs/website: poissonconsulting/poissontemplate
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 70859e68..7aaabbe1 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -1,5 +1,8 @@
-destination: docs
 url: https://poissonconsulting.github.io/extras/
+
+template:
+  package: poissontemplate
+
 reference:
 - title: Manipulate Objects
   desc: Functions to numericise and fill objects
@@ -187,5 +190,4 @@ reference:
   - '`pextreme`'
   - '`sextreme`'
   - '`as_list_unnamed`'
-template:
-  package: poissontemplate
+

From 0ffe00265878af6068ce7cc824dddc48570b2977 Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 19:32:02 -0800
Subject: [PATCH 19/20] add development mode auto for dev website

---
 _pkgdown.yml | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/_pkgdown.yml b/_pkgdown.yml
index 7aaabbe1..572db40d 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -3,6 +3,9 @@ url: https://poissonconsulting.github.io/extras/
 template:
   package: poissontemplate
 
+development:
+  mode: auto
+
 reference:
 - title: Manipulate Objects
   desc: Functions to numericise and fill objects

From d00b9ac12cdff1f4871f1b58ff0afab939dec4fd Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Mon, 6 Jan 2025 19:35:24 -0800
Subject: [PATCH 20/20] move and update installation instructions

---
 README.Rmd | 46 ++++++++++++++++++++++++++++++++--------------
 README.md  | 50 ++++++++++++++++++++++++++++++++++++++------------
 2 files changed, 70 insertions(+), 26 deletions(-)

diff --git a/README.Rmd b/README.Rmd
index dfa46668..2dd2cd7c 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -13,6 +13,7 @@ knitr::opts_chunk$set(
 )
 ```
 
+
 # extras <a href="https://poissonconsulting.github.io/extras/"><img src="man/figures/logo.png" align="right" height="138" alt="extras website" /></a>
 
 <!-- badges: start -->
@@ -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
-
-<!-- To install the latest release from [CRAN](https://cran.r-project.org) -->
-```{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
@@ -107,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 <https://poissonconsulting.github.io/extras/>.
+
+### 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.
@@ -121,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 cb20a035..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 latest release from [CRAN](https://cran.r-project.org) -->
-
-To install the developmental version from
-[GitHub](https://github.com/poissonconsulting/extras)
-
-``` r
-# install.packages("pak")
-pak::pak("poissonconsulting/extras")
-```
-
 ## Demonstration
 
 ### Summarise MCMC Samples
@@ -116,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
+<https://poissonconsulting.github.io/extras/>.
+
+### 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