R/deprecated/acceptance_prob_slow.R

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

#This script calculates the acceptance probability for the MH algorithm for sampling from
#the posterior parameter distribution.

#' @export
#' @import compiler

acceptance_prob <- function(s,         #fixed permutation
                            theta_0,   #current theta
                            theta_1,   #proposed theta
                            prior,     #function to calculate the prior density for any theta
                            net,              #observed graph
                            formula_rhs,    #rhs formula of model
                            ordering = NULL, #ordering,
                            change_on = NULL, #Change on statistics - can be specified if already calculated to speed up
                            change_off =NULL, #Change off statistics
                            ...
){
  if(prior(theta_1) ==0){return(0)}
  if(prior(theta_0) ==0){return(1)}
  
  #define network size
  n <- network.size(net)
  e_obs <- length(net$mel)
  net <- as.BinaryNet(net)
  if(net$isDirected()){e<-n*(n-1)
  }else{e <- (n)*(n-1)/2}
  
  if(length(s)!= e){stop("permutation is the wrong length")}
  
  #make edge list that we can permute
  edges <- combn(seq(1,n),2)
  if(net$isDirected()){
    edges <- as.list(cbind(as.data.frame(edges),as.data.frame(edges[c(2,1),])))
  }
  else{
    edges <- as.list(as.data.frame(edges))
  }
  
  if(is.null(change_on) | is.null(change_off)){
    tmp <- lolog_change_stats(net,s,formula_rhs)
    change_on  <- tmp$change_on
    change_off <- tmp$change_off
  }
  
  change_on_theta_0 = sapply(change_on,function(x){exp(sum(theta_0*x))})
  change_off_theta_0 = sapply(change_off,function(x){exp(sum(theta_0*x))})
  change_on_theta_1 = sapply(change_on,function(x){exp(sum(theta_1*x))})
  change_off_theta_1 = sapply(change_off,function(x){exp(sum(theta_1*x))})
  
  z_t_theta_0 <- change_on_theta_0 + change_off_theta_0
  z_t_theta_1 <- change_on_theta_1 + change_off_theta_1
  
  edges_0 <- edges[s]
  tmp <- sapply(edges_0,function(x){net$getDyads(x[1],x[2])})
  c_t_theta_0 <-  change_on_theta_0*tmp +  change_off_theta_0*(1-tmp)
  
  #commented out since we are not changing permutation
  ##edges_1 <- edges[s]
  ##tmp <- sapply(edges_1,function(x){net$getDyads(x[1],x[2])})
  c_t_theta_1 <-  change_on_theta_1*tmp +  change_off_theta_1*(1-tmp)
  

  return(prod((c_t_theta_1/z_t_theta_1) *(z_t_theta_0/c_t_theta_0)) * (prior(theta_1)/prior(theta_0)))
}

acceptance_prob = cmpfun(acceptance_prob)

#makes pretty much no difference
###### SPEED TESTS ######
#Test

# library(microbenchmark)
# microbenchmark(
#   tmp = acceptance_prob(s_0 = sample(seq(1,120),120,replace = FALSE),
#                         s_1 = sample(seq(1,120),120,replace = FALSE),
#                         theta = model$theta,
#                         net = net_business,
#                         formula_rhs = "edges + triangles + star(c(2,3))")
#   ,times = 1000)