R/deprecated/HMC_Blolog.R

In duncan-clark/Blolog: Bayesian Inference for LOLOG Network Models

#This script carries out HMC for Bayesian LOLOG:

HMC_Blolog <- function(formula_rhs,
                                net,
                                n_perms,
                                verbose = TRUE,
                                prior = function(theta){(1/(10**length(theta)))*(sum(theta>10)==0)*(sum(theta < -10)==0)},
                                burnin = 10,  #Number of burnin iterations for each permutation
                                samples = 10, #Number of samples from each chain
                                samples_post = 10,  #Number of samples to be take from each permutation
                                start = NULL, #specify the start if want the same start for all chains
                                cores =NULL, #if is non null specifies a cluster for parallelisation
                                model = NULL,
                                L_steps = 10,
                                epsilon = 0.1,
                                momentum_sigma = NULL
){
  #define network size
  n <- network.size(net)
  e_obs <- length(net$mel)
  if(get.network.attribute(net,"directed")){e<-n*(n-1)
  }else{e <- (n)*(n-1)/2}
  
  #define momentum_sigma if null
  if(is.null(momentum_sigma)){momentum_sigma <- diag(1,length(start))}
  
  #Get the variational estimates as a starting point for the sampling
  if(is.null(start)){start <- lologVariational(as.formula(paste("net ~",formula_rhs,sep = "")))$theta}
  
  
  posterior_samples <- list()
  length(posterior_samples) <- samples_post*n_perms
  
  sample_func <- function(start){
    s <- sample(1:e,e,replace = FALSE)
    probs <- c()
    
    #Calculate the change stats to speed things up:
    tmp <- lolog_change_stats(net,s,formula_rhs)
    change_on <- tmp$change_on
    change_off <- tmp$change_off
    
    if(verbose){print("Beginning burn in, printing acceptance probabilties")}
    burnin_results <- list()
    burnin_accept <- runif(burnin,0,1)
    if(is.null(start)){theta_old <- start}
    else{theta_old = start}
    if(burnin !=0){
      for(i in 1:burnin){
        
        #allow for specified proposal distribution
        theta_prop <- HMC_proposal(theta_old, 
                                   net,
                                   s,
                                   formula_rhs,
                                   prior = prior,
                                   L_steps = L_steps,
                                   epsilon = epsilon,
                                   momentum_sigma = momentum_sigma,
                                   change_off = change_off,
                                   change_on = change_on)
        
        prob <- acceptance_prob(s,
                                theta_old,
                                theta_prop,
                                prior,
                                net,
                                formula_rhs,
                                change_off = change_off,
                                change_on = change_on)
        
        probs[i] <- prob
        
        if(verbose){print(i);print(prob)}
        accept <- (burnin_accept[i] <= prob)
        if(accept){
          burnin_results[[i]] <- theta_prop
          theta_old <- theta_prop}
        else{burnin_results[[i]] <- theta_old}
      }
    }
    
    #carry out MCMC
    #need to store results
    results <- list()
    results[[1]] <- theta_old
    length(results) <- samples
    if(verbose){print("Beginning sampling, printing acceptance probabilties")}
    samples_accept <- runif(samples,0,1)
    
    for(i in 2:samples){
      
      #allow for specified proposal distribution
      theta_prop <- HMC_proposal(theta_old, 
                                 net,
                                 s,
                                 formula_rhs,
                                 prior = prior,
                                 L_steps = L_steps,
                                 epsilon = epsilon,
                                 momentum_sigma = momentum_sigma,
                                 change_off = change_off,
                                 change_on = change_on)
      
      prob <- acceptance_prob(s,
                              theta_old,
                              theta_prop,
                              prior,
                              net,
                              formula_rhs,
                              change_off = change_off,
                              change_on = change_on)
      
      if(verbose){print(i);print(prob)}
      
      accept <- (samples_accept[i] <= prob)
      if(accept){
        results[[i]] <- theta_prop
        theta_old <- theta_prop}
      else{results[[i]] <- theta_old}
    }
    return(list(results = results[sample(1:samples,samples_post)],
                probs = probs))
  }
  
  if(!is.null(cores)){
    cluster <- makeCluster(cores)
    clusterEvalQ(cl=cluster,library("lolog"))
    clusterEvalQ(cl=cluster,library("statnet"))
    clusterEvalQ(cl=cluster,library("MASS"))
    clusterExport(cl = cluster, varlist = c("acceptance_prob","lolog_change_stats"),envir = globalenv())
    clusterExport(cl = cluster, varlist = c("HMC_proposal","start","net","formula_rhs","e","prior","model","sample_func"),envir = environment()) 
    doParallel::registerDoParallel(cluster)
    posterior_samples <- foreach(j=1:n_perms)%dopar%{sample_func(start)$results}
    stopCluster(cl=cluster)
  }
  else{
    posterior_samples <- foreach(j=1:n_perms)%do%{
      tmp <- sample_func(start)
      if(verbose){
        print("The median acceptance prob was")
        print(median(min(tmp$probs,1)))
      }
      tmp$results}
  }
  posterior_samples <- do.call(c,posterior_samples)
  return(posterior_samples)
}