#### Testing Bayesian functions.
test_that("Functions for Bayesian mixAR work", {
model <- new("MixARGaussian", prob = c(.7,.3), scale = c(1, 2), arcoef = list(-0.5, 1.1))
model2 <- new("MixARGaussian", prob = c(.5,.3, .2), scale = c(1, 2, 3),
arcoef = list(c(-0.5,0.1), 1.1,0.5))
nsim = 20
g <- length(model@prob)
p <- max(model@order)
set.seed(15042020)
y <- mixAR_sim(model, 300, init = 0)
############## Test sampling function for mixing weights
bk1 <- sapply(model@arcoef@a, function(x) 1 - sum(x))
bk2 <- sapply(model2@arcoef@a, function(x) 1 - sum(x))
mu1 <- model@shift / bk1 ; mu2 <- model2@shift / bk2
set.seed(123)
Zpi <- sampZpi(y, model@order, model@prob, mu1,
model@arcoef@a, model@scale, nsim, d = 1)
## Now change the length of "d" and have identical result
set.seed(123)
Zpi2 <- sampZpi(y, model@order, model@prob,mu2,
model@arcoef@a, model@scale, nsim, d = c(1, 1))
expect_identical(Zpi, Zpi2)
#expect_equal_to_reference(Zpi, "sampZpi.rds")
## Expect error if length of "d" is different from 1 or g
expect_error(sampZpi(y, model@order, model@prob, mu1,
model@arcoef@a, model@scale, nsim, d = rep(1, 3)))
expect_error(sampZpi(y, model2@order, model2@prob, mu2,
model2@arcoef@a, model2@scale, nsim, d = rep(1, 2)))
## test dimensions of output
expect_equal(length(Zpi), 3) ## output should be a list of 3
expect_equal(dim(Zpi$mix_weights), c(nsim, g))
expect_equal(dim(Zpi$latentZ), c(length(y) - p, g))
expect_equal(length(Zpi$nk), g)
############## Test sampling function for component mean/shift
#set.seed(234)
mu_shift <- sampMuShift(y, model@order, 1/model@scale^2, Zpi$nk,
model@shift, Zpi$latentZ, model@arcoef@a, nsim = 20)
#expect_equal_to_reference(mu_shift, "sampMuShift.rds")
## test dimension of output
expect_equal(length(mu_shift), 2) ## a list of 2
expect_equal(dim(mu_shift$shift), c(nsim, g))
expect_identical(dim(mu_shift$shift), dim(mu_shift$mu))
############## Test sampling function for component scale/precision
#set.seed(345)
# a = 0.2 , c = 2 are hyperparameters
sigma_tau <- sampSigmaTau(y, model@order, 1/model@scale^2, Zpi$nk, model@arcoef@a,
mu_shift$mu[nsim, ], Zpi$latentZ, a = 0.2 , c = 2, nsim)
#expect_equal_to_reference(sigma_tau, "sampSigmaTau.rds")
## test dimension of output
expect_equal(length(sigma_tau), 3) ## a list of 2
expect_equal(dim(sigma_tau$scale), c(nsim, g))
expect_identical(dim(sigma_tau$scale), dim(sigma_tau$precision))
expect_equal(length(sigma_tau$lambda), nsim)
############## Test bayes_mixAR
nsim = 50
burnin = 20
tau <- c(.1, .2)
#set.seed(456)
bayes <- bayes_mixAR(y, model, fix_shift = FALSE, a = 0.2, c = 2, tau, nsim = 50, burnin = 20)
#expect_equal_to_reference(bayes, "bayes_mixAR.rds")
## test dimension of output
expect_equal(length(bayes), 11) ## a list of 11
expect_equal(dim(bayes$mix_weights), c(nsim - burnin, g))
expect_identical(dim(bayes$mix_weights), dim(bayes$scale))
expect_identical(dim(bayes$mix_weights), dim(bayes$precision))
expect_identical(dim(bayes$mix_weights), dim(bayes$mu))
expect_identical(dim(bayes$mix_weights), dim(bayes$shift))
expect_equal(dim(bayes$LatentZ), c(length(y) - p, g))
expect_equal(length(bayes$ARcoeff), g)
expect_equal(dim(bayes$ARcoeff$Component_1), c(nsim - burnin, model@order[1]))
expect_equal(dim(bayes$ARcoeff$Component_2), c(nsim-burnin, model@order[2]))
expect_equal(bayes$n_samp, nsim-burnin)
expect_equal(bayes$ncomp, g)
expect_error(bayes_mixAR(y, model, fix_shift = FALSE, a = 0.2, c = 2,
tau=c(1, 2, 3), nsim = 11, burnin = 1))
expect_error(bayes_mixAR(y, model2, fix_shift = FALSE, a = 0.2, c = 2,
tau=c(1, 2), nsim = 11, burnin = 1))
############## Test Choose_pk
ar_orders <- Choose_pk(y, model, fix_shift = FALSE, tau, pmax = 5, method = "NULL",
par = NULL, nsim = 10)
expect_equal(ncol(ar_orders$x), g + 1)
expect_equal(sum(ar_orders$x[ , g + 1]), 1)
expect_error(Choose_pk(y, model, fix_shift = FALSE, tau, pmax = 5,
method = c("Poisson", "NULL"), par = NULL, nsim = 10))
expect_error(Choose_pk(y, model, fix_shift = FALSE, tau = c(1,2,3), pmax = 5,
method = "NULL", par = NULL, nsim = 10))
expect_error(Choose_pk(y, model2, fix_shift = FALSE, tau = c(1,2), pmax = 5,
method = "NULL", par = NULL, nsim = 10))
expect_error(Choose_pk(y, model2, fix_shift = FALSE, tau = c(1,2), pmax = 1,
method = "NULL", par = NULL, nsimv= 10))
expect_error(Choose_pk(y, model2, fix_shift = FALSE, tau = c(1,2), pmax = -1,
method = "NULL", par = NULL, nsim = 10))
############## Test marg_loglik
marginal <- marg_loglik(y, model, tau, nsim = 10, prob_mod = ar_orders$x[1, g + 1])
expect_equal(length(marginal), 5)
expect_equal(length(marginal$phi_hd), g)
expect_equal(length(marginal$phi_hd[[1]]), model@order[1])
expect_equal(length(marginal$phi_hd[[2]]), model@order[2])
expect_equal(length(marginal$prec_hd), g)
expect_identical(length(marginal$prec_hd), length(marginal$mu_hd))
expect_identical(length(marginal$prec_hd), length(marginal$weig_hd))
expect_error(marg_loglik(y, model, tau, nsim = 10, prob_mod = -1))
expect_error(marg_loglik(y, model, tau, nsim = 10, prob_mod = 2))
expect_error(marg_loglik(y, model, tau = c(1, 2, 3), nsim = 10,
prob_mod = ar_orders$x[1, g + 1]))
expect_error(marg_loglik(y, model2, tau = c(1, 2), nsim = 10,
prob_mod = ar_orders$x[1, g + 1]))
marg_loglik(y, model2, tau=c(.1), nsim = 10, prob_mod = 0.3)
########### Test label_switch()
label_switch(rbind(bayes$scale[1:25,], bayes$scale[26:30, 2:1]) , m=5)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.