tests/testthat/test-utils.R

test_that("Expansion is correct length.", {
  J <- 5 # Number of basis functions.
  basis <- legendre_basis(5)
  expansion <- expand(1,basis)
  expect_equal(length(expansion),J+1)
})

test_that("Is model matrix correct dimension.", {
  J <- 5; n <- 10
  basis <- legendre_basis(J)
  x <- rep(0,n)
  phi <- modelmat(x,basis)
  expect_equal(dim(phi),c(n,J+1))
})

test_that("Subsets are split evenly?", {
  x <- stats::runif(10,-1,1)
  partition <- split_subset(x,c(-1,1))
  subset1_len <- length(x[ x <= partition[1,2] ])
  expect_equal(subset1_len,5)
  x_odd <- stats::runif(9,-1,1)
  expect_error(split_subset(x_odd,c(-1,1)),
               regexp = "Cannot split an odd number of observations.")
})

test_that("inside() returns correct elements.",{
  x <- 1:10
  x_reduced <- inside(x,c(1,5), end = TRUE)
  expect_equal(x_reduced,1:5)
})

test_that("Correct number of subsets in partition", {
  x <- runif(2^6,-1,1)
  part <- partition(x,3) # Divide subsets in half 3 times
  k <- nrow(part)
  expect_equal(k,2^3)
  #Check that we have the same number of obs in each subset.
  for (i in 1:k){
    if (i == k) x1 <- inside(x,part[i,],end = T) else x1 <- inside(x,part[i,],end = F) # obs inside the i'th subset
    expect_equal(length(x1),2^3)
  }
})

test_that("Partition only works on [-1,1].",{
  x <- runif(2^6,-2,2)
  expect_error(partition(x,3),
               regexp = "Observations must be between -1,1")
})

test_that("Allocate function allocates correctly.", {
  x <- c(5,15,25)
  part <- rbind(
    c(0,10),c(10,20),c(20,30)
  )
  x_alloc <- allocate(x,part)
  expect_equal(x_alloc,c(1,2,3))
})

test_that("Remap keeps obs within [-1,1].", {
  x <- runif(10,-0.5,0.5)
  x1 <- remap(x,c(-0.5,0.5))
  x_outside <- c(x1[x1 > 1], x1[x1 < -1])
  expect_true(length(x_outside) == 0)
})

## WE CREATE SOME SYNTHETIC DATA FOR TEST PURPOSES.
J <- 5
prior <- new_prior(J)
set.seed(100)
beta_true <- MASS::mvrnorm(1,prior@beta,prior@Sigma) # Implies sd = 1
basis <- legendre_basis(J)
fun <- function(x){
  out <- 0.0
  for(j in 1:(J+1)) out <- out + beta_true[j]*basis[[j]](x)
  return(out)
}
n <- 5000
set.seed(100)
x <- runif(n,-1,1)
r <- rep(NA,n)
set.seed(100)
for(t in 1:n) r[t] <- stats::rnorm(1,fun(x[t]),1)

test_that("Posterior converges to true values.", {
  posterior <- update_dist(x,r,prior,basis)
  expect_true(all.equal(posterior@beta,beta_true,tolerance = 0.1))

  # Check convergence of variance parameter.
  var_estimate <- posterior@b / (posterior@a - 1)
  expect_true( all.equal(var_estimate,1,tolerance = 0.1) )
})

test_that("Simulating from correct distribution.", {
  J <- 5
  distribution <- new_prior(J)
  beta_trials <- c(); var_trials <- c()
  set.seed(100)
  for (t in 1:10000){
    sim <- sim(distribution,return_var = TRUE)
    beta_trials <- rbind(beta_trials,sim$beta)
    var_trials <- c(var_trials,sim$var)
  }
  beta_mean <- colMeans(beta_trials)
  var_mean <- mean(var_trials)

  # First test mean and covariance of beta
  expect_true( all.equal(distribution@beta,beta_mean, tolerance = 0.1) )
  exp_var <- distribution@b / (distribution@a - 1)
  expect_true( all.equal(exp_var*distribution@Sigma,cov(beta_trials), tolerance = 0.1) )

  # Test mean and variance of sigma^2
  expect_true(all.equal(exp_var,var_mean, tol = 0.1))
  var_var <- distribution@b^2 / ( (distribution@a - 1)^2 * (distribution@a - 2) ) # True variance of sigma^2
  expect_true(all.equal(var_var,var(var_trials), tolerance = 0.1))
})
dfcorbin/MABsim documentation built on April 26, 2020, 8:26 a.m.