mcse_update_eigen.R

library(cmdstanr)
library(data.table)
library(ggplot2)
library(ggridges)
library(rstan)
library(posterior)
set_cmdstan_path("~/your_cmdstan_path")
# statcomp/benchmarks is working directory
gp_lst = rstan::read_rdump("./gp_pois_regr/gp_pois_regr.data.R")

#base
gp_mod_old = cmdstan_model("./gp_pois_regr/gp_pois_regr.stan")
old_fit = gp_mod_old$sample(
  data = gp_lst, 
  seed = 123, 
  num_chains = 4, 
  num_cores = 4
)

old_fit$cmdstan_summary()

old_stanfit = rstan::read_stan_csv(old_fit$output_files())
print(old_stanfit)


# Switch math via git checkout dependency/eigen-3.3.7 in cmdstan
# and delete old model and .o etc
gp_mod = cmdstan_model("./gp_pois_regr/gp_pois_regr.stan")
fit = gp_mod$sample(
  data = gp_lst, 
  seed = 123, 
  num_chains = 4, 
  num_cores = 4
)

fit$cmdstan_summary()

stanfit = rstan::read_stan_csv(fit$output_files())
foo = print(stanfit, se = TRUE)
old_stan_fit_dt = as.data.table(as.data.frame(old_stanfit))
old_stan_fit_dt[, version := "old"]
stan_fit_dt = as.data.table(as.data.frame(stanfit))
stan_fit_dt[, version := "new"]

full_stan_fit_dt = rbind(old_stan_fit_dt, stan_fit_dt)

mlt_fit_dt = melt(full_stan_fit_dt, id.vars = "version")
sub_mlt_dt = mlt_fit_dt[grepl("f_tilde", variable) & grepl("[5-8]", variable)]
ggplot(sub_mlt_dt,
  aes(x = value, y = variable, color = version)) +
  stat_binline(bins = 50, alpha = 0.25) +
  theme_bw() +
  ggtitle("Histogram of")


ggplot(sub_mlt_dt,
  aes(x = value, y = variable, color = version)) +
  stat_density_ridges(quantile_lines = TRUE, quantiles = c(0.025, 0.5, 0.975), alpha = 0.2) +
  theme_bw() 


fit_ess = stan_ess(stanfit)
old_fit_ess = stan_ess(old_stanfit)

make_mcse_dt = function(fit_mod) {
  fit_names = names(fit_mod)
  par_list = list()[seq_len(length(fit_names))]
  names(par_list) = fit_names
  for (i in seq_len(length(fit_names))) {
    x_par = extract_variable_matrix(fit_mod, fit_names[i])
    mcse_q = mcse_quantile(x_par)
    mcse_dt = data.table(par = fit_names[i], 
      x_mcse_mean = mcse_mean(x_par),
      x_mcse_med = mcse_median(x_par),
      x_mcse_quan_low = mcse_q[1],
      x_mcse_quan_high = mcse_q[2],
      x_mcse_sd = mcse_sd(x_par)
    )
    par_list[[fit_names[i]]] = mcse_dt  
  }
  return(rbindlist(par_list))
}

mcse_fit_dt = make_mcse_dt(stanfit)
mcse_old_fit_dt = make_mcse_dt(old_stanfit)
mcse_diff_dt = as.data.table(as.matrix(mcse_fit_dt[, c(2:6)]) - as.matrix(mcse_old_fit_dt[, c(2:6)]))
mcse_diff_dt[, par := mcse_fit_dt[, 1]]
setcolorder(mcse_diff_dt, "par")
knitr::kable(mcse_diff_dt)