## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----echo = FALSE, eval = TRUE------------------------------------------------
set.seed(123)

## ----setup, eval = TRUE, echo = FALSE-----------------------------------------
library(nrMCmix)

## ----load, echo = TRUE, eval = TRUE-------------------------------------------
data("galaxy", package = "bmixture")

## ----galaxy, echo = TRUE, eval = TRUE-----------------------------------------
model <- specify_model("Gaussian",
                       hyperpar_alpha = TRUE,
                       hyperpar_delta = TRUE,
                       hyperpar_baseline = TRUE)

sampler <- specify_sampler(SAM = TRUE)

fit <- nrMCmix(y = galaxy, model = model, sampler = sampler,
               warm_up = 1000, thin = 10)


## ----print, echo = TRUE, eval = TRUE------------------------------------------
fit
print(summary(fit))

## ----plot, echo = TRUE, eval = TRUE, fig.align='center',fig.width=5, fig.height=4----
plot(fit, n_breaks = 33)

## ----log_p, echo = TRUE, eval = TRUE, fig.align='center',fig.width=5, fig.height=4----
plot(fit,4,par_name = "log_p")

## ----K, fig.align='center',fig.width=5, fig.height=4--------------------------
plot(fit,4, par_name = "K")

## ----densities, fig.align='center',fig.width=5, fig.height=3------------------
par(mfrow = c(1,2))
plot(fit,5,par_name = "alpha", lab_name = "alpha")
plot(fit,5,par_name = "delta", lab_name = "delta")
par(mfrow = c(1,1))

## ----iris---------------------------------------------------------------------
data(iris)

y <- as.matrix(iris[1:4])
fit <- nrMCmix(y, sampler = specify_sampler(SAM = TRUE),
               model = specify_model("MGaussian", dense = FALSE),
               warm_up = 1000, thin = 10, refresh_console = 0)

## ----iris_plot, fig.align='center',fig.width=5, fig.height=5------------------
{
  op <- par(no.readonly = TRUE)
par(mfrow = c(4,4), mar = c(1.1,1.1,1.1,1.1))
for(i in 1:4){
  for(j in 1:4){
    if(i == j){
      plot(fit,idx1 = i, verbose = FALSE, ggplot2 = FALSE)
    }else{
      plot(fit,idx1 = i, idx2 = j, verbose = FALSE, ggplot2 = FALSE)
    }
  }
}
par(op)
}

## ----iris_plot2, fig.align='center',fig.width=5, fig.height=5-----------------
plot(fit,3)

## ----binomial-----------------------------------------------------------------
model <- specify_model("Binomial", alpha0 = 2.1, beta0 = 1.4)

set.seed(123)
y <- generate_data(n = 1000, K = 7, model = model)$y

## ----PY_binomial--------------------------------------------------------------
fit <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(SAM = TRUE),
               warm_up = 1000,
               refresh_console = 0)

## ----binomial_plot, fig.align='center',fig.width=5, fig.height=5--------------
plot(fit)

## ----binomial_comparison------------------------------------------------------
#set the seed
set.seed(123)

#marginal gibbs
MG <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(gibbs = TRUE),
               warm_up = 100, N = 1000,refresh_console = 0)

#non-reversible
NR <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(),
               warm_up = 100, N = 1000,refresh_console = 0)

#non-reversible + split and merge
SAM_NR <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(SAM = TRUE),
               warm_up = 100, N = 1000,refresh_console = 0)

#non-reversible informed + split and merge
SAM_NR_inf10 <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(SAM = TRUE,
                                         informed = TRUE,
                                         m = 10),
               warm_up = 100, N = 1000,refresh_console = 0)

#get the entropy traces
entropy <- cbind(MG = extract_params(MG,"entropy"),
                 NR = extract_params(NR,"entropy"),
                 SAM = extract_params(SAM_NR,"entropy"),
                 SAM_NR_inf10 = extract_params(SAM_NR_inf10,
                                             "entropy"))
      
#get the overall computational cost 
#(number of predictive density evaluations)
cost <- cbind(MG = extract_params(MG,"n_pred_calls"),
              NR = extract_params(NR,"n_pred_calls"),
              SAM = extract_params(SAM_NR,"n_pred_calls"),
              SAM_NR_inf10 = extract_params(SAM_NR_inf10,
                                            "n_pred_calls"))
             

df <- data.frame(log(apply(entropy,2,coda::effectiveSize)),
                 log(apply(cost,2,sum)),
                 log(apply(entropy,2,coda::effectiveSize))-
                   log(apply(cost,2,sum)))
colnames(df) <- c("log(ESS)","log(n° target evaluations)","log(ESS) per target evaluation")
rownames(df) <- c("Marginal Gibbs", "Non-reversible",
                  "Non-reversible + Split and Merge",
                  "Non-reversible-informed (m = 20) + Split and Merge")
knitr::kable(df)