R/deprecated/Blolog.R

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

#This script carries out Bayesian inference on the parameters of the LOLOG model

#' @export

Blolog <- function(formula_rhs,
                   net,
                   n_perms,
                   proposal_step = function(theta,y){MASS::mvrnorm(1,theta,Sigma = diag(1,length(theta)))},
                   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
                   ){
  #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}
  
  #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)
    
    #Calculat 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 <- proposal_step(theta_old,y)
        prob <- acceptance_prob(s,
                                theta_old,
                                theta_prop,
                                prior,
                                net,
                                formula_rhs,
                                change_on=change_on,
                                change_off=change_off)
        
        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 <- proposal_step(theta_old,y)
      theta_prop <- proposal_step(theta_old,y)
      prob <- acceptance_prob(s,
                              theta_old,
                              theta_prop,
                              prior,
                              net,
                              formula_rhs,
                              change_on=change_on,
                              change_off=change_off)
      
      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(results[sample(1:samples,samples_post)])
  }
  
  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("proposal_step","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)}
    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(probs,1)))
      }
      tmp$results}
  }
  posterior_samples <- do.call(c,posterior_samples)
  return(posterior_samples)
}