prob = function(theta,mu,sig2){
  dnorm(log(theta/(1-theta)),mu,sqrt(sig2))/theta/(1-theta)
}

thetas = seq(0.001,0.999,by=0.001)

par(mfrow=c(1,2))

plot(thetas,prob(thetas,0,0.1),
type="l",ylim=c(0,5),xlab=expression(theta),ylab="Density")
lines(thetas,prob(thetas,0,1),type="l",col=2)
lines(thetas,prob(thetas,0,3),type="l",col=3)
lines(thetas,prob(thetas,0,4),type="l",col=4)
lines(thetas,prob(thetas,0,6),type="l",col=6)
legend("topleft",legend=c("sig2=0.1","sig2=1","sig2=3","sig2=4","sig2=6"),col=c(1:4,6),lty=1)
title("mu=0")

plot(thetas,prob(thetas,1,0.1),
type="l",ylim=c(0,10),xlab=expression(theta),ylab="Density")
lines(thetas,prob(thetas,1,1),type="l",col=2)
lines(thetas,prob(thetas,1,3),type="l",col=3)
lines(thetas,prob(thetas,1,4),type="l",col=4)
lines(thetas,prob(thetas,1,6),type="l",col=6)
legend("topleft",legend=c("sig2=0.1","sig2=1","sig2=3","sig2=4","sig2=6"),col=c(1:4,6),lty=1)
title("mu=1")


# Posterior mean (prior 3)
sum(thetas^6*(1-thetas)^3*prob(thetas,0,3)*(thetas[2]-thetas[1]))/
sum(thetas^5*(1-thetas)^3*prob(thetas,0,3)*(thetas[2]-thetas[1]))

# Posterior density (prior 3)
post = thetas^6*(1-thetas)^3*prob(thetas,0,3)
post = post/sum(post*(thetas[2]-thetas[1]))


# Priors and posteriors
par(mfrow=c(1,1))
plot(thetas,prob(thetas,0,3),
type="l",ylim=c(0,3.5),xlab=expression(theta),ylab="Density",col=3)
lines(thetas,dbeta(thetas,1,1))
lines(thetas,dbeta(thetas,4,2),col=2)
lines(thetas,dbeta(thetas,7,5),lwd=3)
lines(thetas,dbeta(thetas,10,6),col=2,lwd=3)
lines(thetas,post,col=3,lwd=3)
points(6/10,0,col=4,pch=16)