#data
x <- 4
n <- 10

#prior
alpha <- 1
beta <- 1

#proposal scale
scale <- 0.005

#start value:
t.t <- 0.9

#Output:
T <- NULL

for(i in 1:1000) {
	#propose new value
	s <- min(t.t * (1 - t.t) / (1 + scale), scale) 

	a <- t.t * (t.t * (1 - t.t) / s - 1)
	b <- (1 - t.t) * (t.t * (1 - t.t) / s - 1)

	tilde.t <- rbeta(1, a, b)	

	#likelihood:
	pi <- dbinom(x, n, tilde.t) / dbinom(x, n, t.t)
	#prior:
	pi <- pi * dbeta(tilde.t, alpha, beta) / dbeta(t.t, alpha, beta)

	#proposal:
	ss <- min(tilde.t * (1 - tilde.t) / (1 + scale), scale) 

	aa <- tilde.t * (tilde.t * (1 - tilde.t) / ss - 1)
	bb <- (1 - tilde.t) * (tilde.t * (1 - tilde.t) / ss - 1)	

	pi <- pi * dbeta(t.t, aa, bb) / dbeta(tilde.t, a, b)

	#make probability
	pi <- min(pi, 1)
	
	if(rbinom(1, 1, pi) == 1) {
		t.t <- tilde.t
	} else {
		t.t <- t.t
	}

	T[i] <- t.t
}

#Inspect the chain
par(mar = c(4,6,1,1), las = 1, cex = 1.7)
plot(T, type ="l", lwd = 2, main = "", xlab = "", ylab = "", axes = FALSE)
axis(1)
axis(2)
mtext(1, line = 2, text = "Iteration / Time", cex = 1.7)
par(las = 0)
mtext(2, line = 4, text = expression(theta), cex = 1.7)

#Correct invariant?
par(mar = c(4,5,1,1), las = 1, cex = 1.7)
hist(T, 100, freq = FALSE, main = "", xlab = "", ylab = "", axes = FALSE)
tt <- seq(min(T), max(T), length.out = 1000)
lines(tt, dbeta(tt, alpha + x, beta + n - x), lwd = 2)
axis(1)
axis(2)
mtext(1, line = 3, text = expression(theta), cex = 1.7)
par(las = 0)
mtext(2, line = 3, text = "Density / Histogram", cex = 1.7)