inst/doc/Tutorial_for_main_functions.R

In XDNUTS: Discontinuous Hamiltonian Monte Carlo with Varying Trajectory Length

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

## ----eval = FALSE-------------------------------------------------------------
# theta0 <- lapply(seq_len(N_chains),myfunction)

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

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

## ----eval = TRUE--------------------------------------------------------------
#observed data
X <- 50

#hyperparameteers
a0 <- 10
b0 <- 10

#list of arguments
arglist <- list(data = X, hyp =  c(a0,b0))

## ----eval = TRUE--------------------------------------------------------------
#function for the negative log posterior and its gradient
#with respect to the continuous components
nlp <- function(par,args,eval_nlp = TRUE){
  
  if(eval_nlp){
    #only the negative log posterior
    
    #overflow
    if(any(par > 30)) return(Inf) 
    
    #conversion of the r value
    r <- ceiling(1 + args$data*plogis(par[2])) - 1
    
    #ensure that r is not zero
    if(r == 0){
      r <- 1
    }
    
    #output
    out <- sum(log(seq_len(r-1))) + 
      (args$data + args$hyp[1] + args$hyp[2])*log(1+exp(-par[1])) + 
      par[1]*(args$data - r + args$hyp[2]) + par[2] + 2*log(1+exp(-par[2]))
    if(r > 2) out <- out - sum(log(seq(args$data - r + 1,args$data - 1)))
    
    return(out)
    
  }else{
    #only the gradient
    
    #overflow
    if(any(par > 30)) return(Inf) 
    
    #conversion of the r value
    r <- ceiling(1 + args$data*plogis(par[2])) - 1
    
    #conversion of the r value
    r <- ceiling(1 + args$data*plogis(par[2])) - 1
    
    #output
    return( (args$data - r + args$hyp[2]) - 
              (args$data + args$hyp[1] + args$hyp[2])*(1-plogis(par[1])) )
  }
  
}

## ----warning=FALSE, message=FALSE, eval = TRUE--------------------------------
#MCMC
set.seed(1)
chains <- xdnuts(theta0 = lapply(1:4,function(x) c(omega = rnorm(1),r_hat = rnorm(1))),
                 nlp = nlp,
                 args = arglist,
                 k = 1,
                 thin = 1,
                 parallel = FALSE,
                 method = "NUTS",
                 hide = TRUE)

## ----eval = TRUE--------------------------------------------------------------
chains

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains, type = 2, gg = FALSE)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains, type = 3)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains, type = 4)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains, type = 5)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains, type = 6)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains, type = 7)

## ----warning=FALSE, message=FALSE, out.width='70%',out.height='70%', eval = TRUE----
summary(chains)

## ----eval = TRUE, dpi = 300---------------------------------------------------
#MCMC
set.seed(1)
chains <- xdnuts(theta0 = lapply(1:4,function(x) rnorm(2)),
                 nlp = nlp,
                 args = arglist,
                 k = 2,
                 thin = 1,
                 parallel = FALSE,
                 method = "XHMC",
                 hide = TRUE,
                 tau = 0.01)

## ----eval = TRUE--------------------------------------------------------------
#define the function to be applied to each sample
f <- function(x,XX) {
  c(
    plogis(x[1]), #inverse logistic for the probability
    ceiling(1 + XX*plogis(x[2])) - 1 #number of trials
  )
}
original_chains <- xdtransform(X = chains, which = NULL,
                               FUN = f,XX = arglist$data,
                               new.names = c("p","r"))

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 2, gg = FALSE)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 4)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 6)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
summary(original_chains)

## ----eval = TRUE--------------------------------------------------------------
data(viscosity)
viscosity

#create the list function
arglist <- list(data =  as.matrix(viscosity[,-1]),
                hyp = c(0, #mean a priori for mu
                        1000, #variance a priori for mu
                        0.5,1,0.5,1 #inverse gamma iperparameters
                        )
                )

## ----eval = TRUE--------------------------------------------------------------
nlp <- function(par,args,eval_nlp = TRUE){
    
    if(eval_nlp){
      #only the negative log posterior

      #overflow
      if(any(abs(par[2:3]) > 30)) return(Inf)
      
      return(par[2] * ( prod(dim(args$data)) + args$hyp[3] ) + 
               (sum( (args$data - par[-(1:3)])^2 ) + 
                  2*args$hyp[4])*exp(-par[2])/2 +
               par[3] * (length(par[-(1:3)]) + 
                           args$hyp[5]) +  
               (sum( (par[-(1:3)] - par[1])^2 ) + 
                  2+args$hyp[6])*exp(-par[3])/2 +
               (par[1] - args$hyp[1])^2/2/args$hyp[2])
      
    }else{
      #only the gradient
      
      #overflow
      if(any(abs(par[2:3]) > 30)) return(rep(Inf,9))
      
      c(
        -sum( par[-(1:3)] - par[1] )*exp(-par[3]) + (
          par[1] - args$hyp[1])/args$hyp[2], #mu derivative
        
        prod(dim(args$data)) + args$hyp[3] - 
          (sum( (args$data - par[-(1:3)])^2 ) +
             2*args$hyp[4])*exp(-par[2])/2, #omega derivative
        
        length(par[-(1:3)]) + args$hyp[5] - 
          (sum( (par[-(1:3)] - par[1])^2 ) + 
             2+args$hyp[6])*exp(-par[3])/2, #omega_a derivative
        
        -apply(args$data - par[-(1:3)],1,sum)*exp(-par[2]) + 
          (par[-(1:3)] - par[1])*exp(-par[3]) #random effects gradient
      )
      
    }
  
}

## ----eval = TRUE--------------------------------------------------------------
#MCMC
set.seed(1)
chains <- xdnuts(theta0 = lapply(1:4,function(x) {
                      out <- rnorm(3 + 6)
                      names(out) <- c("mu","log_s2","log_s2_a",
                                      paste0("mu",1:6))
                      out}),
                 nlp = nlp,
                 args = arglist,
                 k = 0, #no discontinuous components
                 thin = 1,
                 parallel = FALSE,
                 method = "HMC",
                 hide = TRUE,
                 L = 31,
                 control = set_parameters(L_jitter = 5,
                                          N_adapt = 1e3,
                                          keep_warm_up = TRUE))

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(chains,warm_up = TRUE)

## ----eval = TRUE--------------------------------------------------------------
original_chains <- xdtransform(X = chains,which = 2:3,
                               FUN = exp,new.names = c("s2","s2_a"))

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 3)

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 2, which = 1:3, gg = FALSE, scale = 0.5)#fixed

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 2, which = 4:9, gg = FALSE, scale = 0.5)#random

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
plot(original_chains, type = 6)

## ----out.width='70%', out.height='70%', eval = TRUE---------------------------
summary(original_chains, q.val = c(0.05,0.5,0.95))

## ----out.width='70%', out.height='70%', eval = TRUE, dpi = 300----------------
#extract samples into matrix
theta <- xdextract(original_chains,collapse = TRUE)

#compute prediction
y_hat <- sapply(1:6, function(i){
  rnorm(NROW(theta),theta[,3+i],sqrt(theta[,2]))
})

#plot prediction
op <- par(no.readonly = TRUE)
par(mar=c(5.1, 4.1, 2.1, 4.1), xpd=TRUE)
plot(NULL, xlim = c(1,6), ylim = c(15,85), xlab = "Subject",
     ylab =  "Viscosity", bty = "L")
for(i in 1:6){
  #data
  points(rep(i,7),viscosity[i,-1], pch = 16)
  
  #data 95% credible intervals for the prediction
  lines(rep(i,2),quantile(y_hat[,i],c(0.025,0.975)), col = 3, lwd = 3)
  
  #random effects 95% credible intervals
  lines(rep(i,2),quantile(theta[,3+i],c(0.025,0.975)), col = 4, lwd = 4)
}
legend("topright",inset=c(-0.2,-0.2), lty = 1, lwd = 2, col = c(3,4),
       title = "95% Credible Intervals",
       legend = c("prediction","random effects"),
       bty = "n", bg = "transparent", cex = 0.8)
par(op)