R/deprecated/Blolog_QNHMC.R

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

#This script carries out HMC for Bayesian LOLOG:

#' @export

Blolog_QNHMC <- function(formula_rhs,
                       net,
                       n_perms,
                       proposal_step = HMC_proposal,
                       verbose = TRUE,
                       prior = function(theta){(1/(10**length(theta)))*(sum(theta>10)==0)*(sum(theta < -10)==0)},
                       prior_grad = NULL,
                       sample = 10,  #Number of sample iterations for 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
                       diff_perms = TRUE, #if false will use the same perm for each sampling - allows assesment of convergence
                       L_steps = 10,
                       epsilon = 0.1,
                       momentum_sigma = NULL,
                       QNHMC_k = 10, #number of chains in the QNHMC method.
                       ...
){
  t <- proc.time()
  #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}

  s <- sample(1:e,e,replace = FALSE)
  
  #Make the lolog formula
  formula <- as.formula(paste("net ~",formula_rhs,sep = ""))
  
  #Calculate gradient of the log liklihood function for LOLOG, without prior
  log_lik_grad <- function(q,change_on,change_off,statistics = NULL){
    q <- as.vector(q)
    
    
    
    
    
    #derivative of change statistics on top
    if(is.null(statistics)){
      top_deriv <- lolog::calculateStatistics(formula)
    }else{
      top_deriv <- statistics
    }
    
    #derivative of change statistics on bottom
    tmp <- lapply(1:length(change_on),function(i){
      (change_off[[i]]*exp(sum(q*change_off[[i]])) + change_on[[i]]*exp(sum(q*change_on[[i]])))/
        (exp(sum(q*change_off[[i]])) + exp(sum(q*change_on[[i]])))})
    bottom_deriv <- apply(as.matrix(as.data.frame(tmp)),1,sum)
    
    #return their sum with the correct signs
    return(top_deriv  - bottom_deriv)
  }
  
  #gradient of posterior 
  posterior_grad <- function(q,change_on,change_off,statistics = NULL){
    
    if(is.null(prior_grad)){
      prior_deriv <- (numDeriv::grad(prior,q)/prior(q))
    }
    else{prior_deriv <- prior_grad(q)/prior(q)}
    
    if(sum(is.na(prior_deriv)) != 0){
      prior_deriv <- rep(0,length(q))
    }
    #derivative of change statistics on top
    if(is.null(statistics)){top_deriv <- lolog::calculateStatistics(formula)}else{top_deriv = statistics}
    #derivative of change statistics on bottom
    tmp <- lapply(1:length(change_on),function(i){
      (change_off[[i]]*exp(sum(q*change_off[[i]])) + change_on[[i]]*exp(sum(q*change_on[[i]])))/
        (exp(sum(q*change_off[[i]])) + exp(sum(q*change_on[[i]])))})
    bottom_deriv <- apply(as.matrix(as.data.frame(tmp)),1,sum)
    
    #return their sum with the correct signs
    return(-prior_deriv + top_deriv  - bottom_deriv)
  }
  
  #For QNHMC define method to estimate the hessian given a number of previous vectors
  QNHMC_hessian_calc <- function(obs,
                                 theta,
                                 s,
                                 net,
                                 formula_rhs,
                                 change_on,
                                 change_off,
                                 statistics){
    n = length(obs)
    dim = length(obs[[1]])
    grads = lapply(obs,function(x){posterior_grad(x,change_on,change_off,statistics)})
    #Rank the points in terms of there liklihood.
    #Do this by using the acceptance prob function
    #First randomly choose base obs
    base_obs = obs[[sample(1:length(obs),1)]]
    rank <- sapply(obs,function(x){
      acceptance_prob(s,
                      theta_0 = base_obs,
                      theta_1 = x,
                      prior=function(x){1},
                      net = net,
                      formula_rhs,
                      change_off = change_off,
                      change_on = change_on)
    })
    rank <- sort(rank,decreasing = F,index.return = T)$ix
    obs = obs[rank]
    grads = grads[rank]
    
    #need to ensure that h_k is positive deinfite <=> t(s_k))%*%y_k > 0 for all k
    #so remove the x_k value from list if t(s_k))%*%y_k <= 0
    
    #use observed information matrix as hessian as starting point
    #h <- lapply(grads,function(x){outer(x,x)})
    #h <- Reduce("+",h)/length(h)
    tmp <- posterior_grad(theta,change_on,change_off,statistics)  
    h <- outer(tmp,tmp)
    #h <- diag(1,dim)
    for(k in 2:n){
      #print(h)
      if(k<=length(obs)){
        s_k = as.vector(obs[[k]] - obs[[(k-1)]])
        y_k = as.vector(grads[[k]] - grads[[(k-1)]])
        while(t(s_k)%*%y_k >=0){
          print(k)
          print(t(s_k)%*%y_k)
          obs = obs[-k]
          grads = grads[-k]
          s_k = as.vector(obs[[k]] - obs[[(k-1)]])
          y_k = as.vector(grads[[k]] - grads[[(k-1)]])
        }
        rho = (1/t(s_k)%*%y_k)[1,1]
        h = (diag(1,dim) - rho*(y_k%*%t(s_k)) )%*%h%*%(diag(1,dim) - rho*(s_k%*%t(y_k))) + outer(s_k,s_k)
        #h = (diag(1,dim) - rho*(s_k%*%t(y_k)) )%*%h%*%(diag(1,dim) - rho*(y_k%*%t(s_k))) + rho*outer(s_k,s_k)
      }
    }
    
    return(h)
  }
  
  #Need to define different sample func for QNHMC since it uses multiple chains
  sample_func <- function(start){
      
      if(diff_perms){s <- sample(1:e,e,replace = FALSE)}
      
      #Calculate the change stats to speed things up:
      if(verbose){print("Calculating LOLOG change stats")}
      tmp <- lolog_change_stats(net,s,formula_rhs)
      change_on <- tmp$change_on
      change_off <- tmp$change_off
      formula <- as.formula(paste("net ~",formula_rhs,sep = ""))
      statistics <- lolog::calculateStatistics(formula)
      if(verbose){print("Completed calculating LOLOG change stats")}
      
      #Make containers
      sample_results <- lapply(1:sample,function(x){lapply(1:QNHMC_k,function(x){NULL})})
      sample_proposals <- lapply(1:sample,function(x){lapply(1:QNHMC_k,function(x){NULL})})
      sample_probs <- lapply(1:sample,function(x){lapply(1:QNHMC_k,function(x){NULL})})
      sample_accept <- lapply(1:sample,function(x){lapply(1:QNHMC_k,function(x){runif(1,0,1)})})
      #Define the starting thetas
      theta_old <- lapply(1:QNHMC_k,function(x){start()})
      if(sample !=0){
        for(i in 1:sample){
          print(i)
          #need to use for loop since Hessian is iteratively updated:
          for(j in 1:QNHMC_k){
            hessian = QNHMC_hessian_calc(obs = theta_old[-j],
                                         theta = theta_old[[j]],
                                         s=s,
                                         net=net,
                                         formula_rhs=formula_rhs,
                                         change_on = change_on,
                                         change_off = change_off,
                                         statistics = statistics)
            if(verbose){print(paste("The hessian for chain ", j, "was"))
                        print(hessian)}
            tmp = proposal_step(theta_old[[j]], 
                                net,
                                s = s,
                                formula_rhs=formula_rhs,
                                prior = prior,
                                prior_grad = prior_grad,
                                L_steps = L_steps,
                                epsilon = epsilon,
                                momentum_sigma = hessian,
                                change_off = change_off,
                                change_on = change_on)  
            
            sample_proposals[[i]][[j]] = tmp$proposal
            if(verbose){
              print(paste("The proposed step for chain ",j, " was"))
              print(tmp$proposal)
            }
            prob = tmp$prob_factor*acceptance_prob(s,
                                                   theta_old[[j]],
                                                   tmp$proposal,
                                                   prior,
                                                   net,
                                                   formula_rhs,
                                                   change_off = change_off,
                                                   change_on = change_on)
            prob <- as.vector(prob)
            if(verbose){
              print(paste("The step probability for chain ",j, " was"))
              print(prob)
            }
            if(is.nan(prob)){prob <- 0}
            sample_probs[[i]][[j]] <- prob
            if(sample_accept[[i]][[j]] < prob){
              sample_results[[i]][[j]] <- sample_proposals[[i]][[j]]
            }
            else{
              sample_results[[i]][[j]] <- theta_old[[j]]
            }
            theta_old[[j]] <- sample_results[[i]][[j]]
          }
         }
      }
      
    return(list(sample_results = sample_results,
                sample_probs = sample_probs,
                sample_proposals = sample_proposals))
    }

  
  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("log_lik_grad","QNHMC_hessian_calc,","proposal_step","start","net","formula_rhs","e","prior","prior_grad","verbose","sample_func","s","randomwalk_sigma","cov_start"),envir = environment()) 
    doParallel::registerDoParallel(cluster)
    prob_medians <- c()
    posterior_samples <- foreach(j=1:n_perms)%dopar%{
      sample_func(start)}
    stopCluster(cl=cluster)
  }
  else{posterior_samples <- foreach(j=1:n_perms)%do%{sample_func(start)}}
  
  #separate the chains
  sample <- lapply(posterior_samples,function(x){
    lapply(1:QNHMC_k,function(i){
      list(sample_results = lapply(x$sample_results,function(x){x[[i]]}),
           sample_probs = lapply(x$sample_probs,function(x){x[[i]]}),
           sample_proposals = lapply(x$sample_proposals,function(x){x[[i]]}))
      
    })
  })
  
  t <- proc.time() - t
  
  return(list(full_results = sample,time = t[3]))
}