tests/testthat/test-samplers.R

In BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics

context("Test basic functionality of all samplers")

set.seed(1)
library(BayesianTools)

test_that("Samplers work in principle", {
  
  ll = function(x) sum(dnorm(x, log = T))
  setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10))
  
  iter = 200
  iterSMC = 1
  
  samp = getPossibleSamplerTypes()
  
  for(i in 1:(length(samp$BTname))){
    
     if(samp$BTname[i] == "SMC") settings = list(iterations = iterSMC, initialParticles = 100, consoleUpdates = 1e+8) else settings = list(iterations = iter, consoleUpdates = 1e+8)
     invisible(capture.output(suppressMessages(suppressWarnings(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings)))))

      expect_true("bayesianOutput" %in% class(out))
      print(paste("checking" , samp$BTname[i]))
      
    }
    
})



skip_on_cran()
skip_on_ci()

testExact = F

if(testExact){
  iter = 500000
  start = 500
  iterSMC = 400
  library(msm) 
  # library(Matching) # Because this is only needed for manual testing I removed it from
  # the package dependencies. If you are running the manual tests please install
  # the package yourself.
  
  
  test_that("1-d samplers work", {
    
    
    if(test == "exact") skip_on_cran()
    
    ll = function(x) sum(dnorm(x, log = T))
    setup = createBayesianSetup(ll, lower = c(-10), upper = c(10))
    
    samp = getPossibleSamplerTypes()
    
    for(i in 1:length(samp$BTname)){
      
      if(samp$univariatePossible[i] == T){
        if (samp$restartable[i] == T){
          settings = list(iterations = iter/2, adaptive = F, consoleUpdates = 1e+8)
          invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
          invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = out, sampler = samp$BTname[i], settings = settings))))
        }else{
          settings = list(adaptive = F, iterations = iter, consoleUpdates = 1e+8)
          invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
        } 
        
        ks <- ks.boot(getSample(out, numSamples = 10000), rnorm(10000))$ks.boot.pvalue  
        
        # Test that distribution is not significally different from gaussian
        expect_true(ks>0.005)
        
        
      }
      
      
    }
    
    
  })
  
  
  
  test_that("2-D samplers with restarting works", {
    
    ll = function(x) sum(dnorm(x, log = T))
    setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10))
    
    
    samp = getPossibleSamplerTypes()
    
    for(i in 1:(length(samp$BTname)-1)){
      
      if (samp$restartable[i] == T){
        settings = list(iterations = iter/2, consoleUpdates = 1e+8)
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = out, sampler = samp$BTname[i], settings = settings))))
      }else{
        settings = list(iterations = iter, consoleUpdates = 1e+8)
        if(samp$BTname[i] == "SMC") settings = list(iterations = iterSMC, 
                                                    initialParticles = 10000, consoleUpdates = 1e+8)
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
      } 
      
      # Get sample 
      x = getSample(out, numSamples  = 10000)
      y <- rnorm(10000)
      for(z in 1:ncol(x)){
        
        #  ks <- ks.test(x[,z], pnorm)$p.value
        
        ks <- ks.boot(x[,z], y)$ks.boot.pvalue 
        
        # Test that distribution is not significally different from gaussian
        expect_true(ks>0.005)
        
      }
      
    }
    
  })
  
  
  
  test_that("2-d with parallel chains works", {
    
    
    ll = function(x) sum(dnorm(x, log = T))
    setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10), parallel = 2)
    
    
    samp = getPossibleSamplerTypes()
    
    for(i in 1:(length(samp$BTname)-1)){
      
      settings = list(nrChains = 3, iterations = iter, consoleUpdates = 1e+8)
      if(samp$BTname[i] == "SMC") settings = list(nrChains = 3, iterations = iterSMC, consoleUpdates = 1e+8)
      invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings)))) 
      
      
      # Get sample 
      x = getSample(out, numSamples  = 10000)
      y <- rnorm(10000)
      
      for(z in 1:ncol(x)){
        
        #ks <- ks.test(x[,z],pnorm)$p.value  
        ks <- ks.boot(x[,z], y)$ks.boot.pvalue 
        # Test that distribution is not significally different from gaussian
        expect_true(ks>0.005)
        
      }
      
      
    }
    stopParallel(setup)
    
  })
  
  
  test_that("2-d with parallel chains and restart works", {
    
    ll = function(x) sum(dnorm(x, log = T))
    setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10), parallel = 2)
    
    
    samp = getPossibleSamplerTypes()
    
    for(i in 1:(length(samp$BTname)-1)){
      
      if (samp$restartable[i] == T){
        settings = list(nrChains = 3, iterations = iter/2, consoleUpdates = 1e+8)
        if(samp$BTname[i] == "SMC") settings = list(nrChains = 3, iterations = iterSMC, consoleUpdates = 1e+8)
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = out, sampler = samp$BTname[i], settings = settings))))
      }else{
        settings = list(nrChains = 3, iterations = iter, consoleUpdates = 1e+8)
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
      } 
      
      # Get sample and remove first 1000 iterations as burn-in
      x = getSample(out, numSamples  = 10000)
      y <- rnorm(10000)
      
      for(z in 1:ncol(x)){
        
        # ks <- ks.test(x[,z], pnorm)$p.value
        
        ks <- ks.boot(x[,z], y)$ks.boot.pvalue 
        
        # Test that distribution is not significally different from gaussian
        expect_true(ks>0.005)
        
      }
      
    }
    stopParallel(setup)
    
  }
  )
  
  
  
  
  test_that("2-d with external parallelization and restart works", {
    
    # Define function
    FUN <- function(x, sum = TRUE){
      if(sum) sum(dnorm(x, log = T))
      else dnorm(x, log = T)
    }
    
    # Make cluster
    cl <- parallel::makeCluster(2)
    
    # Likelihood function
    ll <- function(pars,...){
      if(is.matrix(pars)) res = parallel::parApply(cl = cl, pars, 1, FUN, ...)
      else res = FUN(pars, ...)
      return(res)
    }
    
    setup <- createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10),
                                 parallel = "external")
    
    samp = getPossibleSamplerTypes()
    
    for(i in 1:(length(samp$BTname)-1)){
      
      if (samp$restartable[i] == T){
        settings = list(nrChains = 1, iterations = iter/2, consoleUpdates = 1e+8)
        if(samp$BTname[i] == "SMC") settings = list(nrChains = 3, iterations = iterSMC, consoleUpdates = 1e+8)
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = out, sampler = samp$BTname[i], settings = settings))))
      }else{
        settings = list(nrChains = 1, iterations = iter, consoleUpdates = 1e+8)
        invisible(capture.output(suppressMessages(out <- runMCMC(bayesianSetup = setup, sampler = samp$BTname[i], settings = settings))))
      } 
      
      # Get sample and remove first 1000 iterations as burn-in
      x = getSample(out, numSamples  = 10000)
      y <- rnorm(10000)
      
      for(z in 1:ncol(x)){
        
        # ks <- ks.test(x[,z], pnorm)$p.value
        
        ks <- ks.boot(x[,z], y)$ks.boot.pvalue 
        
        # Test that distribution is not significally different from gaussian
        expect_true(ks>0.005)
        
      }
      
      
    }
    
    # Stop cluster
    parallel::stopCluster(cl)
    
  }
  )
  
} # if test == "exact"