#install.packages("LaplacesDemon")
#library("LaplacesDemon")

den.t = function(theta,mu,sig,df){
  dt((theta-mu)/sig,df=df)/sig
}
prior.D = function(theta){
  (dnorm(theta,mu.A,sig.A)+
  dnorm(theta,mu.B,sig.B)+
  den.t(theta,mu.C,sig.C,nu))/3
}

sig = 40

mu.A= 900
sig.A = 20

mu.B = 800
sig.B = 80

mu.C = 900
sig.C = sqrt(240)
nu     = 5

pm = (mu.A+mu.B+mu.C)/3


set.seed(125532)
M      = 1000000
thetas = runif(M,600,1000)

par(mfrow=c(1,3))

for (x in c(825,850,875)){

w      = dnorm(thetas,x,sig)*dnorm(thetas,mu.A,sig.A)
ind.A  = sample(1:M,size=M,replace=TRUE,prob=w)
w      = dnorm(thetas,x,sig)*dnorm(thetas,mu.B,sig.B)
ind.B  = sample(1:M,size=M,replace=TRUE,prob=w)
w      = dnorm(thetas,x,sig)*den.t(thetas,mu.C,sig.C,nu)
ind.C  = sample(1:M,size=M,replace=TRUE,prob=w)
w      = dnorm(thetas,x,sig)*(pis[,1]*dnorm(thetas,mu.A,sig.A)+
         pis[,2]*dnorm(thetas,mu.B,sig.B)+pis[,3]*den.t(thetas,mu.C,sig.C,nu))
ind.D  = sample(1:M,size=M,replace=TRUE,prob=w)
a      = c(1,1,1)
pis    = rdirichlet(M,a)
w      = dnorm(thetas,x,sig)*(pis[,1]*dnorm(thetas,mu.A,sig.A)+
         pis[,2]*dnorm(thetas,mu.B,sig.B)+pis[,3]*den.t(thetas,mu.C,sig.C,nu))/ddirichlet(pis,a)
ind.E  = sample(1:M,size=M,replace=TRUE,prob=w)


boxplot(thetas[ind.A],thetas[ind.B],thetas[ind.C],thetas[ind.D],thetas[ind.E],
outline=FALSE,names=c("Post A","Post B","Post C","Post D","Post E"))
abline(h=x,col=2,lwd=3)
points(1:5,c(mu.A,mu.B,mu.C,pm,pm),col=4,pch=16)
title(paste("x=",x,sep=""))
}
legend("bottomright",legend=c("Data point","Prior means"),col=c(2,4),bty="n",lty=2,lwd=3)