R/deprecated/acceptance_prob_slow2.R

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

#Faster acceptance prob:

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) ==0){return(0)}
  if(prior(theta_0) ==0){return(1)}
  
  if(length(s)!= e){stop("permutation is the wrong length")}
  
  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
  }
  
  if(is.null(edge_present)){
    edges <- combn(seq(1,n),2)
    if(get.network.attribute(net,"directed")){
      tmp <- cbind(edges,edges[c(2,1),])
      edges <- lapply(1:length(edges),function(i){tmp[,i]})
    }
    else{
      edges <- lapply(1:length(edges),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])})
  }
  
  change_on_theta_0 = sapply(change_on,function(x){exp(sum(theta_0*x))})
  change_on_theta_1 = sapply(change_on,function(x){exp(sum(theta_1*x))})
  
  
  z_t_theta_0 <- change_on_theta_0 + 1
  z_t_theta_1 <- change_on_theta_1 + 1
  
  c_t_theta_0 <-  change_on_theta_0*edge_present + 1*(1-edge_present)
  c_t_theta_1 <-  change_on_theta_1*edge_present + 1*(1-edge_present)
  
  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)