context("simulate & evaluate")
set.seed(1)
# create fake data
n <- 5E3
nsp <- 3
W <- solve(rWishart(1, 10 + nsp, diag(nsp))[ , , 1])
mu <- matrix(0, n, nsp) # true mean
p = pnorm(mu)
e <- matrix(rnorm(n * nsp), n) %*% chol(W) # true e
z <- mu + e # true z
# Simulate with truncation ------------------------------------------------
# true presence/absence
Y <- ifelse(z > 0, 1, 0)
# Make a bayescomm object with two samples of R and B that match the true
# values
mock_bc = list(
trace = list(
R = matrix(cov2cor(W)[upper.tri(W)], nrow = 2, ncol = 3, byrow = TRUE),
B = replicate(nsp, list(matrix(0, nrow = 2, ncol = 2)))
),
call = list(
X = matrix(1, nrow = n, ncol = 2),
Y = Y
)
)
class(mock_bc) = "bayescomm"
Y_sim = simulate(mock_bc, nsim = 1)[1, , ]
# Empirical probability of a community vector
p_Y = table(as.data.frame(Y)) / nrow(Y)
p_sim = table(as.data.frame(Y_sim)) / nrow(Y)
test_that("simulate.bayescomm works with no environment",{
# How bad an approximation is p_sim for p_Y? Close to zero is good.
kl_divergence = sum(p_Y * (log(p_Y) - log(p_sim)))
expect_less_than(kl_divergence, .001)
})
test_that("likelihood evaluation works with no environment", {
# Empirical log-likelihood from counting communities
ll_empirical = sum(p_Y * log(p_Y))
# log-likelihood from pmvnorm
ll_bc = logLik(mock_bc, mock_bc$call$X[, -1, drop = FALSE], mock_bc$call$Y)
expect_equal(ll_empirical, mean(ll_bc), tol = .01)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.