conjugate distributions

gistfile1.txt

library("R6")
library("purrr")

ConjugateDist <- R6Class("ConjugateDist",
  public = list(
    post_name = NULL,
    post_d = NULL,
    post_r = NULL,
    post_p = NULL,
    post_q = NULL,
    prior_name = NULL,
    prior_d = NULL,
    prior_r = NULL,
    prior_p = NULL,
    prior_q = NULL,
    like_name = NULL,
    like_d = NULL,
    add_prior = function(name, ...) {
      for (i in c("d", "p", "q", "r")) {
        f <- match.fun(paste0(i, name))
        self[[paste0("prior_", i)]] <- partial(f, ..., .lazy = FALSE)
      }
      self$prior_name <- name
    },
    add_post = function(name, ...) {
      for (i in c("d", "p", "q", "r")) {
        f <- match.fun(paste0(i, name))
        self[[paste0("post_", i)]] <- partial(f, ..., .lazy = FALSE)
      }
      self$post_name <- name
    }
  )
)



#' Binomial Likelihood and Beta Prior
#'
#' @param x Number of successes
#' @param size Number of observations
#' @param shape1,shape2 Non-negatative parameters of the Beta distribution
ConjBinomialBeta <- R6Class("BinomialBeta",
  inherit = ConjugateDist,
  public = list(
    initialize = function(x, size, shape1, shape2) {
      shape1_post <- shape1 + sum(x)
      shape2_post <- shape2 + sum(size) - sum(x)
      self$like_d = function(.x) {
        sum(dbinom(x, prob = .x, size = size, log = TRUE))
      }
      self$like_name = "binom"
      self$add_prior("beta", shape1 = shape1, shape2 = shape2)
      self$add_post("beta", shape1 = shape1_post, shape2 = shape2_post)
    })
)

# Beta-Binomial prior
alpha <- 1
beta <- 1
x <- 5
size <- 10
bb <- ConjBinomialBeta$new(x = x, size = size, shape1 = alpha, shape2 = beta)

#
# pois_conj <- function(x, lambda, rate = 1, scale = 1 / rate, ...) {
#   force(rate)
#   n <- length(x)
#   x_sum <- sum(x)
#   size <- k + x_sum
#   prob <- scale / (n * scale + 1)
#   list(
#     dist = "nbinom",
#     params = list(size = size, prob = prob),
#     log_likelihood = function(.x, ...) {
#       sum(dpois(x, lambda = .x, log = TRUE))
#     },
#     prior = dist_funs("gamma", rate = rate),
#     posterior = dist_funs("nbinom", size = size, prob = prob)
#   )
# }
#
# # Poisson likelihood, Gamma Prior with rate or scale
# pois_conj <- function(x, lambda, rate = 1, scale = 1 / rate, ...) {
#   force(rate)
#   n <- length(x)
#   x_sum <- sum(x)
#   shape_post <- k + x_sum
#   scale_post <- scale / (n * scale + 1)
#   list(
#     dist = "gamma",
#     params = list(shape_post = shape_post, scale_post = scale_post),
#     log_likelihood = function(.x) {
#       sum(dpois(x, lambda = .x, log = TRUE))
#     },
#     prior = dist_funs("gamma", shape = shape, rate = rate),
#     posterior = dist_funs("gamma", shape = shape_post, rate = rate_post)
#   )
# }
#
# # normal-gamma
# #\left.\left({\frac {\mu _{0}}{\sigma _{0}^{2}}}+{\frac {\sum _{i=1}^{n}x_{i}}{\sigma ^{2}}}\right)\right/\left({\frac {1}{\sigma _{0}^{2}}}+{\frac {n}{\sigma ^{2}}}\right),
# #{\displaystyle \left({\frac {1}{\sigma _{0}^{2}}}+{\frac {n}{\sigma ^{2}}}\right)^{-1}} \left#({\frac {1}{\sigma _{0}^{2}}}+{\frac {n}{\sigma ^{2}}}\right)^{-1}
#
# # Normal gamma
# norm_mu_norm_conj <- function(x, tau, mu0, tau0, ...) {
#   n <- length(x)
#   mu_post <- (tau0 * mu0 + tau * sum(x)) / (tau0 + n * tau)
#   sigma_post <- 1 / sqrt(tau0 + n * tau)
#   sigma_prior <- 1 / sqrt(tau)
#   list(
#     dist = "norm",
#     param = list(mean = mu_post, sd = sigma_post),
#     log_likelihood = function(.x) {
#       sum(dnorm(x, mean = .x, sd = sigma_prior, log = TRUE))
#     },
#     prior = dist_funs("norm", mean = mu0, sd = sigma_prior),
#     posterior = dist_funs("norm", mean = mu0, sigma_post)
#   )
# }
#
# # Normal. known maen mu.
# norm_tau_gamma_conj <- function(x, mu, shape = 1, rate = 1, scale = 1 / rate, ...) {
#   n <- length(x)
#   shape_post <- shape + n * 0.5
#   rate_post <- rate + (sum(x - mu) ^ 2) * 0.5
#   list(
#     dist = "gamma",
#     param = list(shape = shape_post, rate = rate_post, scale = 1 / rate_post),
#     log_likelihood = function(.x) {
#       sum(dnorm(x, mean = mu, sd = .x, log = TRUE))
#     },
#     prior = dist_funs("gamma", shape = shape, rate = rate),
#     posterior = dist_funs("gamma", shape = shape_post, rate = rate_post)
#   )
# }
#
# # Normal. known maen mu.
# norm_mu_tau_conj <- function(x, shape = 1, rate = 1, scale = 1 / rate, ...) {
#   n <- length(x)
#   shape_post <- shape + n * 0.5
#   rate_post <- rate + (sum(x - mu) ^ 2) * 0.5
#   list(
#     dist = "gamma",
#     param = list(shape = shape_post, rate = rate_post, scale = 1 / rate_post),
#     log_likelihood = function(.x) {
#       sum(dnorm(x, mean = mu, sd = .x, log = TRUE))
#     },
#     prior = dist_funs("gamma", shape = shape, rate = rate),
#     posterior = dist_funs("gamma", shape = shape_post, rate = rate_post)
#   )
# }
#
# # VGAM
# unif_pareto_conj <- function(x, max, scale = 1, shape, ...) {
#   scale_post <- base::max(c(x, max))
#   shape_post = shape + length(x)
#   list(
#     dist = "pareto",
#     params = list(scale = scale_post, shape = shape_post),
#     log_likelihood = function(.x) {
#       sum(dunif(x, min = 0, max = .x, log = TRUE))
#     },
#     prior = dist_funs("pareto", scale = scale, shape = shape),
#     posterior = dist_funs("pareto", scale = scale_post, shape = shape_post)
#   )
# }
#
#
# # Probability grid
# # like ppoints, but optional start min and max
# pgrid <- function(.n = 100, .pmin = NULL, .pmax = NULL) {
#   # values come from ppoints
#   a <- if (n <= 10) 3 / 8 else 1 / 2
#   denom <- (.n + 1 - 2 * a)
#   .pmin <- .pmin %||% 1 / denom
#   .pmax <- .pmax %||% .n / denom
#   seq(.pmin, .pmax, length.out = .n)
# }
#
# qgrid <- function(f, .n = 100, ..., .pmin = NULL, .pmax = NULL) {
#   match.fun(f)(pgrid(.n, .pmin, .pmax), ...)
# }
#
# dgrid <- function(f, .n = 100, ..., .pmin = NULL, .pmax = NULL, qfun = NULL) {
#   qfun <- qfun %||% match.fun(gsub("^d", "q", substitute(f)))
#   x <- qgrid(qfun, .n, ..., .pmin = .pmin, .pmax = .pmax)
#   match.fun(f)(x, ...)
# }
#
# dnormgamma <- function(x, shape = 1, rate = 1, mean = 0, sd = 1, log = FALSE) {
#   tau <- 1 / sd ^ 2
#   p <- (stats::dnorm(x, mean = mean, sd = sd, log = TRUE) +
#         stats::dgamma(tau, shape = shape, rate = rate, log = TRUE))
#   if (!log) p <- exp(x)
#   p
# }
#
# rnormgamma <- function(n, shape = 1, rate = 1, mean = 0, sd = 1, log = FALSE) {
#   tau <- stats::rgamma(n, shape = shape, rate = rate)
#   mu <- rnorm(n, mean = mean, sd = 1 / sqrt(tau))
#   cbind(mu, sd)
# }

Some scratch code for Bayesian conjugate distributions. If I do this, I should probably build off of the distr packages, which already have objects with proporties of most distributions.