cjbayesian / mhadaptive Goto Github PK

set.seed(666^3)
M <- 1E4
K <- 31
Mu <- rnorm(K, mean = 0, sd = 10)
P <- qr.Q(qr(matrix(rnorm(K^2), K)))
SS <- crossprod(P, P*(K:1))

target <- function(x){
  mvtnorm::dmvnorm(x, mean = Mu, sigma = SS, log = TRUE)
}
########
indep.samples <- MASS::mvrnorm(n = M, mu = Mu, Sigma = SS)


MH1 <- mcmc::metrop(target, initial = rnorm(K), nbatch = M)
MH1.samples <- MH1$batch

MH2 <- MHadaptive::Metro_Hastings(li_func = target, 
                                  pars = rnorm(K),
                                  burn_in = 1000,
                                  iterations = M + round(.1 * M))
MH2
MH2.samples <- MH2$trace

library(rstan)
mvmodel <- '
data{
  int <lower=0> k;
  matrix[k , k] Sigma;
  vector[k] mu;
}
parameters{
  vector[k] x;
}
model{
 x ~ multi_normal(mu, Sigma);
}
'
stan.data <- list(
  k = K,
  mu = Mu,
  Sigma = SS
)
StanRun <- stan(model_code = mvmodel, data = stan.data, iter = M, chains = 2 )
Stan.samples <- extract(StanRun, 'x')$x

Means <- list(
  indep = data.frame(est = colMeans(indep.samples), truth = Mu, sampler = "indep"),
  MH = data.frame(est = colMeans(MH1.samples), truth = Mu, sampler = "MH"),
  adapMH = data.frame(est = colMeans(MH2.samples), truth = Mu, sampler = "adaptMH"),
  Stan = data.frame(est = colMeans(Stan.samples), truth = Mu, sampler = "Stan")
)

####
## mean squared error and mean relative error
lapply(Means, function(dt) mean((dt$truth-dt$est)^2))
lapply(Means, function(dt) mean((dt$truth-dt$est)/dt$truth))

###
mean.results <- do.call(rbind, Means)
mean.results$sampler <- factor(mean.results$sampler, 
                               levels = sort(levels(as.factor(mean.results$sampler)), FALSE)  )

library(ggplot2)
p0 <- ggplot(mean.results,
             aes(x  =  sampler, y = (truth-est)/truth,
                 fill = sampler, colour = sampler) ) + 
  stat_summary(fun.data = mean_se, size = 1.5) +
  scale_x_discrete("Sampler") +
  scale_y_continuous(expression((mu-hat(mu))/mu),
                     expand = c(0, 0)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() + 
  NULL
p0

p1 <- ggplot(mean.results,
             aes(x  =  sampler, y = (truth-est)/truth,
                 fill = sampler) ) + 
  geom_boxplot(position = "dodge") +
  scale_x_discrete("Sampler") +
  scale_y_continuous(expression((mu-hat(mu))/mu),
                expand = c(0, 0)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() + 
  NULL
p1

p2 <- ggplot(mean.results,
             aes(x  =  truth, y = est,
                 fill = sampler, colour = sampler) ) + 
  geom_point() +
  scale_x_continuous(expression(mu)) +
  scale_y_continuous(expression(hat(mu))) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  theme_bw() + 
  NULL
p2