# com24a of October 2018,
# doing Nils Exercise 22 b for stk 4021,
# when teaching for Celine 

nn <- 25
mu0 <- 2.345
sigma0 <- 1.234

yy <- mu0 + sigma0*rnorm(nn) # these are the data
ybar <- mean(yy)
Q0 <- sum((yy-ybar)^2) 

# since the idea is to fix a given dataset: 
ybar <- 2.3139
Q0 <- 37.5011 

logL <- function(para)
{
mu <- para[1]
sigma <- para[2]
# 
-nn*log(sigma) - 0.5*(Q0 + nn*(mu-ybar)^2)/sigma^2 
}

minuslogL <- function(para)
{-logL(para)}

nils <- nlm(minuslogL,c(0,1),hessian=T)
ML <- nils$estimate
Jhat <- nils$hessian
se <- sqrt(diag(solve(Jhat)))
showme <- cbind(ML,se)
print(round(showme,4))

# ML     se 
# 2.3139 0.2450 mu 
# 1.2248 0.1733 sigma 

big <- 135

logprior <- function(para)
{
mu <- para[1]
sigma <- para[2]
# 
ok <- 1*(abs(mu) <= 5.00)*(abs(log(sigma)) <= 10.00)
1*ok*(1/sigma) - big*(1 - ok)
}

burnin <- 1000 
sim <- burnin + 10^4 # with 1000 as run-in 
mcmc <- 0*(1:sim) %*% t(1:2)  

mcmc[1, ] <- c(4,3)
# mcmc[1, ] <- c(ML[1],ML[2])

uu <- runif(sim)
ok <- 0*(1:sim)
pracc <- 0*(1:sim) 

for (ss in 2:sim)
{
old <- mcmc[ss-1, ]
prop <- old + 0.22*se*rnorm(2) 
# 
aux11 <- logprior(prop) - logprior(old) 
aux12 <- logL(prop) - logL(old)  
# 
oi <- exp(aux11+aux12)
pracc[ss] <- min(c(1,oi))
ok[ss] <- 1*(uu[ss] <= pracc[ss])
mcmc[ss, ] <- (1-ok[ss])*old + ok[ss]*prop
} 

plot(mcmc) # all, including burn-in time

mcmcstar <- mcmc[(burnin+1):sim, ] # where I omit the first burnin period

par(bty="l")
par(cex.lab=1.4142) 

plot(mcmcstar,
   xlab="simulated mu",
   ylab="simulated sigma") 

## the output can be seen as realisations from
#  pi(mu, sigma | data) 
#  and can be used to address all questions
#  regarding (mu, sigma)