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)