# Edpsy 590BAY
# Fall 2019
# c.j.anderson
#
# Function beta_binomial:   Impact of prior on posterior
#
#  Uses a grid method
#
#  Note 
#
# Input y = number of "sucesses"
#       n = number of trials
#       shape.1 = first shape parameter of beta
#       shape.2 = second shaper parameter of beta
#
#  Note: 
#   o beta(1,1) is the uniform distribution
#   o shape.1 < shape.2 is right skewed
#   o shape.1 > shape.2 is left skewed
#   o shape.1 = shape.2 is uni-modal & symetric around 0.5
#   o shape.1 = shape.2 and as shape.1 & shape.2 get larger, then smaller std
#

beta.binomial <- function(y,n,shape.1,shape.2){
	N<- 10000                               # how many p's to check
	p <-  seq(from=0,to=1,length.out=N)     # possible values for probability of sucess
	likelihood <- dbinom(y,n,prob=p)        # binomial likelihood

	layout(matrix(c(1,1,0,0,1,1,3,3,2,2,3,3,2,2,0,0),4,4,byrow=FALSE))

	# Beta prior
	prior.beta <- dbeta(p,shape.1,shape.2)              # beta(shape.1,shape.2) prior
	posterior.beta <- prior.beta*likelihood             # posterior w/ beta(shape.1,shape.2) prior
	posterior.beta <- posterior.beta/sum(posterior.beta)

	plot(p,prior.beta, type="l", main=paste("Beta(",shape.1,",",shape.2,") Prior")     )
	plot(p,likelihood, type="l", main=paste("Binomial Likelihood, y=", y, "& n=", n)   )
	plot(p,posterior.beta, type="l", main=paste("Posterior with Beta(",shape.1,",",shape.2,") Prior") )
}

#
# or All in one plot
#
library(LearnBayes)

BetaBinomial.tri <- function(y,n,shape.1,shape.2){
	triplot(c(shape.1,shape.2),c(y,(n-y)),where="topright")
}