################################################################
# NORMAL MODEL AND NORMAL PRIOR => NORMAL POSTERIOR
################################################################

# Physicist A
m.A  = 900
sd.A = 20

# Physicits B
m.B  = 800
sd.B = 80

# Model & data
sd.X = 40
x    = 850

# Posterior Physicist A
pi.A   = sd.X^2/(sd.X^2+sd.A^2)
mean.A = pi.A*m.A + (1-pi.A)*x
var.A  = pi.A*sd.A^2

# Posterior Physicist B
pi.B   = sd.X^2/(sd.X^2+sd.B^2)
mean.B = pi.B*m.B + (1-pi.B)*x
var.B  = pi.B*sd.B^2


curve(dnorm(x, mean=m.A,sd=sd.A), n=1000,from=500,to=1100, col="blue",
xlab="Measurement",ylab="Density")
curve(dnorm(x, mean=m.B,sd=sd.B), n=1000,from=500,to=1100, col="red",add=TRUE)
points(x,0,pch=16,cex=2)
curve(dnorm(x, mean=mean.A,sd=sqrt(var.A)), n=1000,from=500,to=1100, col="blue",add=TRUE,lty=2)
curve(dnorm(x, mean=mean.B,sd=sqrt(var.B)), n=1000,from=500,to=1100, col="red",add=TRUE,lty=2)

legend("topleft",legend=c("Prior A", "Prior B","Posterior A","Posterior B"),
col=c("blue","red","blue","red"),lty=c(1,1,2,2),bty="n")

legend("topright",legend=c("Data point"),pch=16,bty="n")



################################################################
# POISSON MODEL AND GAMMA PRIOR => GAMMA POSTERIOR
################################################################

# Informative prior 
E  = 1.5      # mean
S2 = 1.0^2    # variance
a  = E^2/S2
b  = E/S2
M  = (a-1)/b  # Mode


# Non-informative prior 
a1 = 0.0001
b1 = 0.0001

par(mfrow=c(1,1))
curve(dgamma(x,a,b), n=1000,from=0,to=6, col="blue",
xlab="Minutes",ylab="Density",ylim=c(0,1.25))

curve(dgamma(x,a1,b1), n=1000,from=0,to=6, col="red",add=TRUE)


# The data
y = c(0,0,0,1,2,3)
n = length(y)

# Posteriors
a2 = a + sum(y)
b2 = b + n
E2 = a2/b2
M2 = (a2-1)/b2

a3 = a1 + sum(y)
b3 = b1 + n
E3 = a3/b3
M3 = (a3-1)/b3

points(y,rep(0,n),pch=16)

curve(dgamma(x,a2,b2), n=1000,from=0,to=6,col="blue",add=TRUE,lty=2)
curve(dgamma(x,a3,b3), n=1000,from=0,to=6,col="red",add=TRUE,lty=2)

legend("topright",legend=c("Informative prior","Uninformative prior"),col=c("blue","red"),lty=1,bty="n")


par(mfrow=c(2,3),mai=c(0.5,0.5,0.2,0.2))
for (i in 1:n){
  curve(dgamma(x,a,b),n=1000,from=0,to=6,xlab="",ylab="",ylim=c(0,1.5),lwd=2)
  a2 = a + sum(y[1:i])
  b2 = b + length(y[1:i])
  curve(dgamma(x,a2,b2),n=1000,from=0,to=6,col=i+1,add=TRUE,lwd=2)
  points(y[1:i],rep(0,i),pch=16,col=2:(i+1))
  mtext(side=1,line=2,"Minutes")
  mtext(side=2,line=2,"Density")
  legend("topright",legend=c(
  paste("Prior: Gamma(",a,",",b,")",sep=""),paste("Posterior after ",i," obs.",sep="")),
  col=c(1,i+1),lty=1,bty="n",lwd=2)
}