R/beta_arm.R

In pikachuuu108/StatComp18018: 18018's final project for 'Statistical Computing' course

#' @title A beta distribution sampler using R
#' @description a function to generate a random sample of size n from the Beta(a, b) distribution by the acceptance-rejection method
#' @param n the number of samples
#' @param a the first parameter of function 'beta'
#' @param b the second parameter of function 'beta'
#' @return a random sample of size \code{n}
#' @import ggplot2
#' @examples
#' sample_beta <- beta_arm(1000,3,2)
#' head(sample_beta,20)
#' sample_beta_theo <- rbeta(1000,3,2)
#' data_beta <- data.frame(sample = c(sample_beta,sample_beta_theo),class = rep(c('empirical','theoretical'),each = 1000))
#' ggplot(data_beta,aes(x = sample,fill = class))+geom_histogram(position="identity", alpha=0.5)
#' @export
beta_arm <- function(n,a,b){
  beta_pdf <- function(x){
    (1/beta(a,b))*(x)^(a-1)*(1-x)^(b-1)#define the pdf of beta distribution
  }
  j <- 0#the time of iteration
  k <- 0#the available sample
  y <- numeric(n)#a zero vector of length n
  while (k < n){
    u <- runif(1)
    j <- j + 1
    x <- runif(1) #random variate from g
    if ((a-1)*x^(a-2) * (b-1)*(1-x)^(b-2) > u) {
      #accept x
      k <- k + 1
      y[k] <- x
    }
  }
  return(y)
}