inst/tree_sampling/tree-sample-test3.R

In skgallagher/TBornotTB: Analyze Tuberculosis Patients in Maryland

## SKG
## 2/4/20
## Calculating likelhood  and making trees
##


devtools::load_all()


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)


## Make samples of data
n_trials <- 100
cluster_summary <- bark

n <- cluster_summary$n
n_pos <- cluster_summary$n_pos
n_vec <- n
n_pos_vec <- n_pos
par <- c("beta_0" = 0, "beta_1" = 0)




## Sample trees of n people with n_pos
sampled_trees <- sample_many_groups(n_vec = n_vec,
                                    n_pos_vec = n_pos_vec,
                                    rep_vec = rep(n_trials, length(n_vec)))

sum_trees <- summarize_cond_trees(sampled_trees)

fast_out <- loglike_obs_data(par,
                             data = cluster_summary,
                                         sampled_trees_summary = sum_trees,
                                         return_neg = FALSE)

fast2 <- optim(par = par,
               fn = loglike_obs_data,
               data = cluster_summary,
               sampled_trees_summary = sum_trees,
               return_neg = TRUE, hessian = TRUE)

fast2$par
sqrt(diag(solve(fast2$hessian)))

## Test
##
parts <- partitions::restrictedparts(n = 8, m = 4, include.zero = FALSE)
permuted_parts <- matrix(unlist(apply(parts, 2, unique_perm)),
                         nrow = 4)
unif_p <- 1 / nrow(permuted_parts)
multi_p1 <- dmultinom(x = c(2, 2, 2, 2), prob = rep(1/4, 4))
multi_p2 <- dmultinom(x = c(5, 1, 1, 1), prob = rep(1/4, 4))
c(unif_p, multi_p1, multi_p2)