R/tests/RMHMC_boot_test.R

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

#This script tests that the new RMHMC boot proposal gives the same as the old one:
library(Blolog)
library(statnet)
library(lolog)
library(coda)
library(doParallel)
library(ggplot2)

setwd("C:/Users/Duncan/Documents/Academics/UCLA_Academics/Networks/LOLOG_Cateloging/1_9_UniversityEmails")
load("processed_data.RData")

n <- get.network.attribute(net,"n")
e <- choose(n,2)

formula_rhs <- "edges+triangles+star(c(2,3))"
prior = function(theta){1}
prior_grad = function(theta){0}
q = c(0,0,0,0)
L_steps = 1
epsilon =1
resamples = 10

RMHMC_boot_1 <- RMHMC_boot_proposal
RMHMC_boot_2 <- RMHMC_boot_proposal_fast
  
RMHMC_boot_test <- function(RMHMC_boot_1,RMHMC_boot_2,times = 10){
  
  s <- sample(1:e,e)
  change_on = Blolog::lolog_change_stats(net,s,formula_rhs)
  change_off = lapply(1:length(change_on),function(x){rep(x,length(change_on[[1]]))})
  
  t_1 <- microbenchmark(RMHMC_boot_1_result <- RMHMC_boot_1(theta_0 = q,
                                                              net = net,
                                                              s = s, 
                                                              formula_rhs = formula_rhs,
                                                              prior = prior,
                                                              prior_grad = prior_grad, 
                                                              L_steps  =1,
                                                              epsilon = 1,
                                                              change_on =change_on,
                                                              resamples = resamples),times = times)
  t_2 <- microbenchmark(RMHMC_boot_2_result <- RMHMC_boot_2(theta_0 = q,
                                                              net = net,
                                                              s = s, 
                                                              formula_rhs = formula_rhs,
                                                              prior = prior,
                                                              prior_grad = prior_grad, 
                                                              L_steps  =1,
                                                              epsilon = 1,
                                                              change_on =change_on,
                                                              resamples = resamples),times = times)
  
  print(paste("Method 2 is ",summary(t_1)$mean/summary(t_2)$mean, " times faster than method 1",sep=""))
  return()
}

RMHMC_boot_test(RMHMC_boot_proposal,RMHMC_boot_proposal_fast,times = 20)