# Bayesian inference problem
# likelihood x|theta ~ N(theta,1)
# prior theta ~ Cauchy
# we observe x = 0.8
# the posterior distribution is intractable
# Bayesian estimate, the posterior expectation is also intractable
# we use the rejection sampling scheme to simulate from the posterior
# we approximate the posterior expectation using the empirical mean over the simulations
# that is a Monte Carlo method

x <- 0.8
target.tilde <- function(theta) exp(x*theta-theta^2/2)/(1+theta^2)
proposal.tilde <- function(theta) 1/(1+theta^2) # Cauchy (Student(1))
M.tilde <- exp(x^2/2)
n <- 10000
simu <- rep(0,n)
nsimu <- 0
for (i in 1:n)
{
bool <- TRUE
while (bool)
{
y <- rt(1,1) # one simulation from a Cauchy
u <- runif(1) # one simulation from a uniform on [0,1]
rho <- target.tilde(y)/(proposal.tilde(y)*M.tilde) # acceptance probability
if (u <= rho) {simu[i] <- y ; bool <- FALSE}
nsimu <- nsimu +1 
}
}