# Edps 590BAY
# Fall 2019
# C.J. Anderson
#
# Figures showing prior x likelihood -> posterior
#

#
# Set up
#
# Need from values for proportion/probability
p <- seq(0,1,length.out=1000)

# To put plots where I want them
layout(matrix(c(1,1,0,0,1,1,3,3,2,2,3,3,2,2,0,0),4,4,byrow=FALSE))

#
# Examples of posteriors for different priors and likelihoods
#

# no impact of uniform prior
#
prior <- dbeta(p,1,1)
plot(p,prior,type='l',main="Prior is beta(1,1)=uniform" )

likelihood <- dnorm(p,.5,.2)
plot(p,likelihood,type='l',main='Likelihod is N(.5,.2)')
posterior <- prior*likelihood/sum(prior*likelihood)

plot(p,posterior,type='l', main='Posterior')


#
# step function for prior
#
zeros <- rep(0,200)
ones <- rep(1,800)
prior <- c(zeros,ones) 

plot(p,prior,type='l',main="Prior is Step Function" )

likelihood <- dnorm(p,.5,.2)
plot(p,likelihood,type='l',main='Likelihod is N(.5,.2)')
posterior <- prior*likelihood/sum(prior*likelihood)

plot(p,posterior,type='l', main='Posterior')

#
# 2 steps 
#
ones <- rep(.4,200)
two <- rep(1,600)
prior <- c(zeros,ones,two)
plot(p,prior,type='l',main="Prior is Step Function" )

likelihood <- dnorm(p,.5,.2)
plot(p,likelihood,type='l',main='Likelihod is N(.5,.2)')
posterior <- prior*likelihood/sum(prior*likelihood)

plot(p,posterior,type='l', main='Posterior')


#
# nice average of distributions to get posterior
#
prior <- dbeta(p,5,2)
plot(p,prior,type='l',main="Prior is beta(5,2)" )

likelihood <- dnorm(p,.5,.2)
plot(p,likelihood,type='l',main='Likelihod is N(.5,.2)')
posterior <- prior*likelihood/sum(prior*likelihood)

plot(p,posterior,type='l', main='Posterior')


#
# average of distributions to get posterior but small variance
#
prior <- dbeta(p,5,2)
plot(p,prior,type='l',main="Prior is beta(5,2)" )

likelihood <- dnorm(p,.5,.01)
plot(p,likelihood,type='l',main='Likelihod is N(.5,.01)')
posterior <- prior*likelihood/sum(prior*likelihood)

plot(p,posterior,type='l', main='Posterior')


#
# average of distributions to get posterior but small variance
#  --- prior and likelihood in opposite directions
#  NOT in lecture notes
#
prior <- dbeta(p,2,5)
plot(p,prior,type='l',main="Prior is beta(5,2)" )

likelihood <- dnorm(p,.5,.2)
plot(p,likelihood,type='l',main='Likelihod is N(.5,.2)')
posterior <- prior*likelihood/sum(prior*likelihood)

plot(p,posterior,type='l', main='Posterior',xlab="Parameter or Theta")