# Chapter 4 R examples

library(bayesrules)
library(tidyverse)
library(janitor)

# Run multiplot function from the Chapter 3 R examples first:


#####
## Run this multiplot function code ONCE:
# Multiplot function for later: just run this once:
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
## End multiplot function code
#############


# Bechtel test example:

# Plots of different beta prior distributions:

b5.11 <- plot_beta(5,11)+ggtitle("Beta(5,11)")
b1.1 <- plot_beta(1,1)+ggtitle("Beta(1,1)")
b14.1 <- plot_beta(14,1)+ggtitle("Beta(14,1)")

# I like the 3x1 array of plots a little better here:
multiplot(b5.11,b1.1,b14.1, cols=1)
#multiplot(b5.11,b1.1,b14.1, cols=2)
#multiplot(b5.11,b1.1,b14.1, cols=3)

# Import data
data(bechdel, package = "bayesrules")

# Take a sample of 20 movies
set.seed(84735)
bechdel_20 <- bechdel %>% 
  sample_n(20)
  
# the first 3 rows of our data set:
bechdel_20 %>% 
  head(3)

# Table of passes and fails for our sample of 20
  
bechdel_20 %>% 
  tabyl(binary) %>% 
  adorn_totals("row")
  
# Plotting the posteriors corresponding to each prior:

b14.22 <- plot_beta(14,22)+ggtitle("Beta(14,22)")
b10.12 <- plot_beta(10,12)+ggtitle("Beta(10,12)")
b23.12 <- plot_beta(23,12)+ggtitle("Beta(23,12)")

multiplot(b14.22,b10.12,b23.12, cols=1)

# Visualizing how the prior and likelihood combine to create the posterior for each prior choice:

p1<-plot_beta_binomial(alpha = 5, beta = 11, y = 9, n = 20)+ggtitle("Beta(5,11) prior")
p2<-plot_beta_binomial(alpha = 1, beta = 1, y = 9, n = 20)+ggtitle("Beta(1,1) prior")
p3<-plot_beta_binomial(alpha = 14, beta = 1, y = 9, n = 20)+ggtitle("Beta(14,1) prior")

multiplot(p1,p2,p3, cols=1)


# Visualizing how the prior and likelihood combine to create the posterior
# for three different data sets, each with sample proportion around 0.46,
# and each with a Beta(14,1) prior, but with very different sample sizes:

d1<-plot_beta_binomial(alpha = 14, beta = 1, y = 6, n = 13)+ggtitle("Small sample")
d2<-plot_beta_binomial(alpha = 14, beta = 1, y = 29, n = 63)+ggtitle("Large sample")
d3<-plot_beta_binomial(alpha = 14, beta = 1, y = 46, n = 99)+ggtitle("Very Large sample")

multiplot(d1,d2,d3, cols=1)


# Plot of Beta(10,90) prior:

b10.90 <- plot_beta(10,90)+ggtitle("Beta(10,90)"); b10.90

# Pessimistic prior overwhelmed by a LARGE amount of data:

plot_beta_binomial(alpha = 10, beta = 90, y = 900, n = 1000)

# Posterior summary for this last example:
summarize_beta_binomial(alpha = 10, beta = 90, y = 900, n = 1000)