rm(list=ls())

#-----------------------------------------------------#
#---   Classical and Bayesian inference (AMS132)   ---#
#---   Class 4: Prior and Posterior distribution   ---#
#-----------------------------------------------------#

#-- The posterior distribution: Example: Lifetime of components
data=c(0.92, 0.72, 1.27, 0.69, 1.97)
n=length(data)

#   b)   prior and posterior distributions:
a=5 ; b=1; a/b; a/b^2
curve(dgamma(x, shape=a, rate=b), 0,8, ylim=c(0,1.5), lwd=3, col="black", 
      xlab=expression(theta), ylab="", main="prior and posterior distr.") # prior 1
curve(dgamma(x, shape=a+n, rate=b+sum(data)), 0,8, add=TRUE, lwd=3, col="blue") # posterior

a=1 ; b=0.2; a/b; a/b^2
curve(dgamma(x, shape=a, rate=b), 0,8, ylim=c(0,1.5), lwd=3, col="black", 
      xlab=expression(theta), ylab="", main="prior and posterior distr.") # prior 2
curve(dgamma(x, shape=a+n, rate=b+sum(data)), 0,8, add=TRUE, lwd=3, col="green") # posterior

#   priors:
curve(dgamma(x, shape=5, rate=1), 0,8, ylim=c(0,1.5), lwd=3, col="blue", 
      xlab=expression(theta), ylab=expression(xi(theta)), main="prior") # prior 1
curve(dgamma(x, shape=1, rate=0.2), 0,8, ylim=c(0,1.5), lwd=3, col="green", add=TRUE) # prior 2

#   posteriors:
curve(dgamma(x, shape=5+n, rate=1+sum(data)), 0,8, ylim=c(0,1.5), lwd=3, col="blue", 
      xlab=expression(theta), ylab="posterior", main="posterior")# posterior 1
curve(dgamma(x, shape=1+n, rate=0.2+sum(data)), 0,8, add=TRUE , lwd=3, col="green")# posterior 2


#-- Likelihood function: Example: Lifetime of components (cont.)
curve(x^n*exp(-x*sum(data)), 0, 8, xlab=expression(theta), ylab="likelihood", lwd=3)

a=1 ; b=0.2; 
curve(x^n*exp(-x*sum(data))*10, 0,8, xlab=expression(theta), lwd=3, ylim=c(0,1))# likelihood
curve(dgamma(x, shape=a, rate=b), 0,8, ylim=c(0,1.5), lwd=3, col="green", add=TRUE) # prior 
curve(dgamma(x, shape=a+n, rate=b+sum(data)), 0,8, add=TRUE, lwd=3, col="blue") # posterior

a=5 ; b=1; 
curve(x^n*exp(-x*sum(data))*10, 0,8, xlab=expression(theta), lwd=3, ylim=c(0,1))# likelihood
curve(dgamma(x, shape=a, rate=b), 0,8, ylim=c(0,1.5), lwd=3, col="green", add=TRUE) # prior 
curve(dgamma(x, shape=a+n, rate=b+sum(data)), 0,8, add=TRUE, lwd=3, col="blue") # posterior