## Okapi R demonstration ## AG Schissler ## Simple Bayes rule calculation ## Posterior grid approximation for binomial experiment with 9 trials and 6 "successes" start_time <- proc.time() ## 1. Define a univariate grid size <- 9 x <- 6 grid_length <- 2e8 p_grid <- seq( from = 0, to = 1, length.out = grid_length) ## 2. Define a flat prior (uniform) prior <- rep(1, grid_length) ## 3. Compute likelihood at each value of the grid likelihood <- dbinom(x = x, size = size, prob = p_grid) ## 4. Compute the product of likelihood and prior unnormalized_posterior <- likelihood * prior ## 5. Normalize (standardize) the posterior to sum to 1 posterior <- unnormalized_posterior / sum( unnormalized_posterior ) ## 6. Summaries ## maximum a posteri? ## p_grid[which.max(posterior)] ## Visualize results (if possible) ## plot( p_grid, posterior, type = "b" ) ## sample from the posterior and summarize samples <- sample(p_grid, size = grid_length, replace = T, prob = posterior) summary(samples) sd(samples) ## compare to mle (p_hat_mle <- x/size) ## sample proportion (mean) (se_mle <- sqrt( size * (p_hat_mle) * (1 - p_hat_mle))) (end_time <- proc.time() - start_time)