R/bernoulli_beta_model.R

In eliobartos/misc: Miscellaneous functions

#' Run bayesian beta - bernoulli model for estimating proportion of 0-1 variable
#'
#' Runs a bayesian estimation of proportion using bernoulli distribution as likelihood and beta distribution as conjugate prior.
#' Posterior distribution is beta distribution.
#'
#' @param alpha Parameter for prior distribution representing the number of success
#' @param beta Parameter for prior distribution representing the number of fails
#' Prior used is Beta(alpha + 1, beta + 1)
#' @param success Number of success cases in your data
#' @param total Total number of cases in your data
#' @param sample_size Size of sample from posterior distribution
#'
#' @return Vector of samples from posterior distribution
#'
#' Posterior distribution is Beta(alpha + 1 + success, beta + 1 + total - success)
#'
#' @examples
#' # No prior information, prior is uniform
#' post = bernoulli_beta_model(0, 0, 20, 100)
#'
#' # Prior succes rate is around 5% with estimation strenght as it was estimated on a sample of 100
#' post2 = bernoulli_beta_model(5, 95, 3, 50)
#'
#' mean(post)
#' quantile(post, probs = c(0.05, 0.95)) # 90% highest density posterior interval
#'
#' mean(post2)
#' quantile(post2, probs = c(0.05, 0.95))
#'
#' @importFrom stats rbeta
#' @export
#' @author Elio Bartoš
bernoulli_beta_model = function(alpha, beta, success, total, sample_size = 1e5) {
  rbeta(sample_size, alpha + 1 + success, beta + 1 + total - success)
}