Nothing
#------------------------------- Beta-binomial -------------------------------#
# Calculate sufficient statistics
binomial_data <- function(data, prior) {
if (!is.matrix(data) & !is.data.frame(data)) {
stop("beta_binom: data must be a matrix or a data frame")
}
if (ncol(data) != 2) {
stop("beta_binom: data must have 2 columns")
}
if (any(!is_wholenumber(data))) {
stop("beta_binom: the data must be whole numbers")
}
n <- data[, 2]
if (any(n <= 0)) {
stop("beta_binom: the values in data column 2 (n) must be positive")
}
y <- data[, 1]
if (any(y < 0)) {
stop("beta_binom: the values in data column 1 (y) must be non-negative")
}
if (any(y > n)) {
stop("beta_binom: in data column 1 (y) cannot exceed column 2 (n), rowwise")
}
# If a default improper prior is used then check for posterior propriety
if (is.character(prior)) {
if (prior == "bda" | prior == "default") {
if (all(pmin(y, n - y) == 0)) {
stop("''bda'': improper posterior unless 0 < y < n at least once")
}
}
}
# Calculate the values needed in beta_binom_marginal_loglik
ny <- n - y
tab_y <- table(y[y > 0])
tab_ny <- table(ny[ny > 0])
tab_n <- table(n)
y_vals <- as.numeric(names(tab_y))
w_y <- as.numeric(tab_y)
y_mat <- cbind(y_vals, w_y)
ny_vals <- as.numeric(names(tab_ny))
w_ny <- as.numeric(tab_ny)
ny_mat <- cbind(ny_vals, w_ny)
n_vals <- as.numeric(names(tab_n))
w_n <- as.numeric(tab_n)
n_mat <- cbind(n_vals, w_n)
#
return(list(y_mat = y_mat, ny_mat = ny_mat, n_mat = n_mat))
}
# Marginal log-likelihood
# posterior for (alpha, beta) not including the prior for (alpha, beta)
# Note: we need to be careful to avoid underflow when either or both
# alpha and beta are very large
beta_binom_marginal_loglik <- function(x, y_mat, ny_mat, n_mat) {
if (any(x <= 0)) {
return(-Inf)
}
s <- x[1] + x[2]
mu <- x[1] / s
f1 <- function(y) {
return(sum(log(mu + (y - 1:y) / s)))
}
t1 <- sum(y_mat[, 2] * vapply(y_mat[, 1], f1, 0))
f2 <- function(ny) {
return(sum(log(1 - mu + (ny - 1:ny) / s)))
}
t2 <- sum(ny_mat[, 2] * vapply(ny_mat[, 1], f2, 0))
f3 <- function(n) {
return(sum(log(1 + (n - 1:n) / s)))
}
t3 <- sum(n_mat[, 2] * vapply(n_mat[, 1], f3, 0))
return(t1 + t2 - t3)
}
# Obvious coding - used only by testthat to test that
# beta_binom_marginal_loglik is correct
check_beta_binom_marginal_loglik <- function(x, y, n) {
if (any(x <= 0)) {
return(-Inf)
}
return(sum(lbeta(y + x[1], n - y + x[2]) - lbeta(x[1], x[2])))
}
# Sample from the conditional posterior distribution of the population
# parameters given the hyperparameters and the data
beta_binom_cond_sim <- function(x, data, n_sim) {
alpha <- x[, 1]
beta <- x[, 2]
y <- data[, 1]
n <- data[, 2]
len_y <- length(y)
theta_sim_vals <- matrix(NA, ncol = len_y, nrow = n_sim)
for (i in 1:len_y) {
theta_sim_vals[, i] <- stats::rbeta(n_sim, alpha + y[i],
beta + n[i] - y[i])
}
colnames(theta_sim_vals) <- paste("p[",1:len_y,"]", sep = "")
return(list(theta_sim_vals = theta_sim_vals))
}
# --------------------------- sim_pred_beta_binom --------------------------- #
#' Simulate from a beta-binomial posterior predictive distribution
#'
#' Simulates \code{nrep} draws from the posterior predictive distribution
#' of the beta-binomial model described in \code{\link{hef}}.
#' This function is called within \code{\link{hef}} when the argument
#' \code{nrep} is supplied.
#' @param theta_sim_vals A numeric matrix with \code{nrow(data)} columns.
#' Each row of \code{theta_sim_vals} contains binomial success probabilities
#' simulated from their posterior distribution.
#' @param data A 2-column numeric matrix: the numbers of successes in column 1
#' and the corresponding numbers of trials in column 2.
#' @param nrep A numeric scalar. The number of replications of the original
#' dataset simulated from the posterior predictive distribution.
#' If \code{nrep} is greater than \code{nrow(theta_sim_vals)} then
#' \code{nrep} is set equal to \code{nrow(theta_sim_vals)}.
#' @return A numeric matrix with \code{nrep} columns. Each column contains
#' a draw from the posterior predictive distribution of the number of
#' successes.
#' @examples
#' rat_res <- hef(model = "beta_binom", data = rat)
#' rat_sim_pred <- sim_pred_beta_binom(rat_res$theta_sim_vals, rat, 50)
#' @export
sim_pred_beta_binom <- function(theta_sim_vals, data, nrep) {
nrep <- min(nrep, nrow(theta_sim_vals))
# Extract the first nrep rows from the posterior sample of thetas
thetas <- theta_sim_vals[1:nrep, , drop = FALSE]
# Function to simulate one set of data
bin_fn <- function(x) {
return(stats::rbinom(n = length(data[, 2]), size = data[, 2], prob = x))
}
return(apply(thetas, 1, bin_fn))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.