tests/testthat/test-weighted-sampling.R

In skgallagher/TBornotTB: Analyze Tuberculosis Patients in Maryland

test_that("likelihood_weighted", {
  par <- c(0, 0)
  prev_par <- c(0,0)
  n <- 5
  x <- 3
  B <- 10
  out <- likelihood_weighted(par = par,
                             prev_par = prev_par,
                             n = n, x = x,
                             B = B)
  expect_true(length(out) == 1)
})


test_that("likelihood_weighted_total", {

  par <- c(0, 0)
  prev_par <- par
  freq <- c(1, 1, 2)
  n_vec <- c(5, 3, 2)
  x_vec <- c(2, 1, 2)
  B <- 100
  out <- loglike_weighted_total(par = par, prev_par = prev_par,
                                freq = freq,
                                n_vec = n_vec, x_vec = x_vec,
                                B = B)
  expect_equal(length(out), 1)
})




# test_that("putting it all together", {
#
#   ## Generate partition list
#   part_list <- generate_part_list(n = 25)
#   g_weight_list <- get_weight_list(part_list)
#
#
#   K <- 10000
#   inf_params <- c("beta_0" = -2, "beta_1" = 1.5)
#   max_size <- 26
#   smear_pos_prob <- .5
#
#   forest <- simulate_flip_til_failure(K = K,
#                                       inf_params = inf_params,
#                                       smear_pos_prob = smear_pos_prob,
#                                       max_size = max_size)
#   bark <- summarize_clusters(forest)
#
#   ################################
#   init_par <- c(-1.82, 1.3)
#   freq <- bark$freq
#   n_vec <- bark$n
#   x_vec <- bark$n_pos
#   B <- 1000
#   iters <- 100
#   out <- iter_est_ws(init_par = init_par,
#                      freq = freq,
#                      n_vec = n_vec,
#                      x_vec = x_vec,
#                      B = B,
#                      iters = iters)
#   out$par
#   sqrt(diag(solve(out$hessian)))
#
#
#   ## Sample trees of n people with n_pos
#   n_trials <- 100
#   n_vec2 <- rep(bark$n, each = n_trials)
#   n_pos_vec <- rep(bark$n_pos, each = n_trials)
#   par <- c("beta_0" = 1, "beta_1" = 0)
#
#
#   sampled_trees2 <- simulate_many_cond_bp(K = length(n_vec2),
#                                           n_vec = n_vec2, n_pos_vec = n_pos_vec,
#                                           part_list = part_list,
#                                           g_weight_list = g_weight_list,
#                                           one_init = TRUE)
#   sum_trees2 <- summarize_cond_trees(sampled_trees2)
#
#
#
#
#   fast2 <- optim(par = par,
#                  fn = loglike_obs_data,
#                  data = bark,
#                  sampled_trees_summary = sum_trees2,
#                  return_neg = TRUE, hessian = TRUE)
#   fast2$par
#   sqrt(diag(solve(fast2$hessian)))
#
#   x <- 7
#   y <- 5
#   out <-  sum_trees2 %>% dplyr::filter(.data$n == x, .data$n_pos == y)
#   hist(out$i_pos, prob = TRUE)
#   my_x <- 0:(x-1)
#   my_y <- dbinom(my_x, size = x-1, p = y / (x-1))
#   lines(my_x, my_y)
#   table(out$i_pos) / sum(table(out$i_pos))
#   my_y
#
#   ## Thoughts
#   ## I think the weights are off
#
# })