# com41a of october 2021
# simple mcmc illustrations 

ff <- function(theta)
{
# trying out different densities, to see how mcmc works: 
# (1/pi) / (1 + theta^2)
# dnorm(theta) 
0.5*dnorm(theta-1.50) + 0.5*dnorm(theta+1.50) 
# exp(-abs(theta)^0.333)
# 1*(abs(theta)>0.10)*1/abs(theta)^3.333
}

curve(ff,-7,7) 
integrate(ff,-100,0)
integrate(ff,0,100)

sim <- 1*10^5
mcmc <- 0*(1:sim)

uu <- runif(sim) # for all the ok = 1 or ok = 0 things
vv <- runif(sim,-1,1) # for the proposal distribution 
pracc <- 0*(1:sim)
ok <- 0*(1:sim) 
delta <- 0.75 # value set after a bit of trial and error 

# i start my mcmc in some position,
# before i let it run its course: 
mcmc[1] <- 7.345 

for (ss in 2:sim) 
{
old <- mcmc[ss-1]
prop <- old + delta*vv[ss]
pracc[ss] <- min(c(1,ff(prop)/ff(old)))
ok[ss] <- 1*(uu[ss] <= pracc[ss])
# 
mcmc[ss] <- ok[ss]*prop + (1-ok[ss])*old
}

plot(mcmc) 

plot(density(mcmc))

# acceptance rate: 
print(c(mean(pracc),mean(ok)))