chilis-1.R

library("tidyverse")
# theta is a parameter representing the probability that a remodeled store sees at least a 5% increase in sales
# Two bosses: Boss1 Theta = 20% and Boss2 Theta = 70%
# we have two competing models (i.e. the Boss1 and Boss2) of the world
# Now make the vector of theta values that represent the models of the world that we have:
thetaVals = c(0.2,0.7)

#CALCULATE PRIOR
# pTheta is the vector of prior probabilities on the theta values.
pTheta = c(0.5,0.5) # Makes a uniform belief distribution.

#SPECIFY OBSERVED DATA
# Specify the data. The follwoing are 3 successes and 9 failure
# and 6 success and 4 failures.
#Data = c(1,1,1,0,0,0,0,0,0,0,0,0)
#Data = c(rep(1,6), rep(0,4))
Data = c(0)  ## one store is a failure
nSuccess = sum( Data == 1 )    ##count # of successful stores
nFail = sum( Data == 0 )    ##count # of failure stores

#CALCULATE LIKELIHOOD
# Compute the likelihood of the data for each value of theta:
pDataGivenTheta = thetaVals^nSuccess * (1-thetaVals)^nFail

#CALCULATE POSTERIOR PROBABILITIES
# Compute the posterior:
pData = sum( pDataGivenTheta * pTheta )
pThetaGivenData = pDataGivenTheta * pTheta / pData   # This is Bayes' rule!

#make data frame for plotting
plotDF = data.frame(thetaVals, prior = pTheta,likelihood = pDataGivenTheta, posterior = pThetaGivenData)

#make tidy data to use facet grid by ProbType for plot
tidyPlotDF = plotDF %>%
  gather("ProbType","Probability",-thetaVals) %>%
  mutate(ProbType = fct_relevel(ProbType, c("prior","likelihood","posterior")))# oreder levels for plot


#create named vector for facet labels
labels = c(prior = "PRIOR\np(theta)",likelihood = "LIKELIHOOD\np(Data|theta)",posterior="POSTERIOR\np(theta|Data)")

#create plot
tidyPlotDF %>%
  ggplot(aes(x = thetaVals, y = Probability, fill = ProbType)) +
  geom_col(width = 0.02, show.legend = FALSE) +
  coord_cartesian(xlim = c(0,1)) + ## create coordinate system
  facet_grid(rows = vars(ProbType),
             labeller = labeller(ProbType = labels),
             scales = "free_y") +
  geom_text(aes(y = Probability * 1.2, label = signif(Probability,2))) +
  labs(title = paste0("Belief Evolution -- ",nFail," failure(s) -- ",nSuccess," success(es)")) +
  theme_minimal(16)