# Sample from a 50/50 mixture of two Gaussians
# at +1/-1, both with variance 1
rmix <- function(n) {
	rnorm(n,2*rbinom(n,1,0.5)-1,1)
}

# Density of the mixture
dmix <- function(x) {
	0.5*dnorm(x,1,1)+0.5*dnorm(x,-1,1)
}

# Give the posterior probability of the +1 mode
postmix <- function(z) {
	1/(1+exp(-2*z))
}

# Posterior predictive density
dpostmix <- function(x,z) {
	p = postmix(z)
	return(p*dnorm(x,1,1)+(1-p)*dnorm(x,-1,1))
}

# Log likelihood ratio of posterior predictive
# to true density
LRpostmix <- function(x,z) {
	log(dmix(x)/dpostmix(x,z))
}

# Relative entropy of posterior predictive to
# true density, by monte carlo
relent.1 <- function(z,m=1e5) {
	mean(LRpostmix(rmix(m),z))
}

# Vectorized version of previous
relent <- function(z,m=1e5) {
	sapply(z,relent.1,m=m)
}

# x = c(0,rmix(1000))
# z = cumsum(x)
# posteriors = postmix(z)
# rels = relent(z)
# par(mfcol=c(3,1))
# plot(0:1000,z,xlab="n",ylab="z",type="l"); abline(h=0,lty=2)
# plot(0:1000,posteriors,xlab="",ylab=expression(P[n](theta==+1)),type="l");abline(h=0.5,lty=2)
# plot(0:1000,rels,xlab="",ylab=expression(D[KL](P[0],P[n])),type="l")