# Data
n = 10
y = c(0,0,0,1,0,1,0,1,0,1)
s = sum(y==1)

# PART I: prior for theta is a mixture of two betas
# -------------------------------------------------

# Prior hyperparameters
a  = 2
b  = 5
c  = 5
d  = 2
pi = 0.5

# B(a,b)
B = function(a,b){
  gamma(a)*gamma(b)/gamma(a+b)
}

# Posterior hyperparameters
a1  = a + s
b1  = b + n - s
c1  = c + s
d1  = d + n -s 
pi1 = 1/(1+(1-pi)/pi*B(a,b)*B(c1,d1)/B(c,d)/B(a1,b1))

theta = seq(0,1,length=1000)

par(mfrow=c(1,2))
plot(theta,dbeta(theta,a,b),xlab=expression(theta),ylab="Prior density",type="l",lwd=3)
lines(theta,dbeta(theta,c,d),col=2,lwd=2)
lines(theta,pi*dbeta(theta,a,b)+(1-pi)*dbeta(theta,c,d),col=3,lwd=2)
legend("top",legend=c("Beta(a,b)","Beta(c,d)","Mixture"),col=1:3,lwd=2,bty="n")
title(paste("(pi,a,b,c,d)=(",pi,",",a,",",b,",",c,",",d,")",sep=""))

plot(theta,dbeta(theta,a1,b1),xlab=expression(theta),
     ylab="Posterior density",type="l",lwd=3)
lines(theta,dbeta(theta,c1,d1),col=2,lwd=2)
lines(theta,pi1*dbeta(theta,a1,b1)+(1-pi1)*dbeta(theta,c1,d1),col=3,lwd=2)
legend("topright",legend=c("Beta(a1,b1)","Beta(c1,d1)","Mixture"),col=1:3,lwd=2,bty="n")
title(paste("(pi1,a1,b1,c1,d1)=(",round(pi1,2),",",a1,",",b1,",",c1,",",d1,")",sep=""))

legend("topleft",legend=c(paste("n=",n,sep=""),paste("s=",s,sep="")),bty="n")


# PART II: Gaussian prior for log odds, i.e. 
# gamma = log(theta/(1-theta))
# ------------------------------------------
M       = 100000
gammas  = rnorm(M,0,1)
thetas  = 1/(1+exp(-gammas))
w       = thetas^s*(1-thetas)^(n-s)
ind     = sample(1:M,replace=TRUE,prob=w)
thetas1 = thetas[ind]
gammas1 = gammas[ind]

par(mfrow=c(2,2))
hist(gammas,prob=TRUE,main="",xlab=expression(gamma))
hist(gammas1,prob=TRUE,main="",xlab=expression(gamma))
hist(thetas,prob=TRUE,main="",xlab=expression(theta),xlim=c(0,1),ylim=c(0,3))
lines(theta,dbeta(theta,a,b),col=2,lwd=2)
hist(thetas1,prob=TRUE,main="",xlab=expression(theta),xlim=c(0,1),ylim=c(0,3.5))
lines(theta,dbeta(theta,a1,b1),col=2,lwd=2)