R/acceptance_prob.R

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

#Faster acceptance prob:

#' @export
#' 
#' 

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
                            edge_present = NULL, #vector giving whether an edge is present or not
                            ...
){
  if(prior(theta_1,net,change_on) ==0){return(0)}
  if(prior(theta_0,net,change_on) ==0){return(1)}
  
  if(is.null(change_on) | is.null(change_off)){
    change_on <- lolog_change_stats(net,s,formula_rhs)
    change_off <- lapply(1:length(change_on),function(x){rep(0,length(change_on[[1]]))})
  }
  
  if(length(s)!= length(change_on)){stop("permutation is the wrong length")}
  
  if(is.null(edge_present)){
    tmp <- combn(seq(1,n),2)
    if(get.network.attribute(net,"directed")){
      tmp <- cbind(tmp,tmp[c(2,1),])
      edges <- lapply(1:dim(tmp)[2],function(i){tmp[,i]})
    }
    else{
      edges <- lapply(1:dim(tmp)[2],function(i){tmp[,i]})
    }
    net_tmp <- as.BinaryNet(net)
    edges <- edges[s]
    edge_present <- sapply(edges,function(x){net_tmp$getDyads(x[1],x[2])})
  }
  
  #tranform edge_present to be -1 when edge present, +1 when no edge 
  edge_present <- 2*(1-edge_present) -1
  
  ratio <- mapply(change_on,edge_present,FUN = function(x,ind){
    (1 + exp(ind*sum(theta_0*x))) / (1 + exp(ind*sum(theta_1*x)))
  })
  
  return(prod(ratio)*(prior(theta_1,net,change_on)/prior(theta_0,net,change_on)))
}