context("updog")
test_that("bin_post works", {
set.seed(18)
ncounts <- c(11, 13)
ssize <- c(21, 37)
prior <- 6
postprob <- bin_post(ncounts = ncounts, ssize = ssize, prior = prior)
expect_equal(sum(postprob), 1)
## test zeros
postprob1 <- bin_post(ncounts = 0, ssize = 2, prior = 2, seq_error = 0)
postprob2 <- bin_post(ncounts = 2, ssize = 2, prior = 2, seq_error = 0)
temp <- ((1 - seq(0, 2) / 2) ^ 2) / 3
expect_true(all(abs(temp / sum(temp) - postprob1) < 10 ^ -14))
expect_true(all(abs(temp / sum(temp) - postprob2[3:1]) < 10 ^ -14))
}
)
test_that("updog works", {
set.seed(734)
ocounts <- c(11, 13, 0, 100)
osize <- c(21, 37, 2, 101)
p1counts <- c(11, 9, 7, 7)
p1size <- c(23, 29, 17, 15)
p2counts <- c(7, 9, 11)
p2size <- c(11, 9, 11)
ploidy <- 6
seq_error <- 0.01
uout <- updog(ocounts = ocounts, osize = osize, p1counts = p1counts,
p1size = p1size, p2counts = p2counts, p2size = p2size,
ploidy = ploidy, seq_error = 0.01, update_od_param = FALSE)
uout2 <- updog(ocounts = ocounts, osize = osize,
ploidy = ploidy, seq_error = 0.01, update_od_param = FALSE)
expect_true(all(abs(rowSums(uout$postmat) - 1) < 10 ^ -14))
# uout <- updog(ocounts = ocounts, osize = osize, p1counts = p1counts,
# p1size = p1size, p2counts = p2counts, p2size = p2size,
# ploidy = ploidy, do_mcmc = TRUE, iterate = TRUE)
#
# plot(uout$opostprob, uout$m_opostprob)
# abline(0, 1)
#
# plot(uout$p1postprob, uout$m_p1postprob)
# abline(0, 1)
#
# plot(uout$p2postprob, uout$m_p2postprob)
}
)
test_that("up_fix increases the up_obj at every iteration", {
dat <- readRDS("dat.RDS")
seq_error <- 0.01
ocounts <- dat[1:(nrow(dat) - 12), 1]
osize <- rowSums(dat[1:(nrow(dat) - 12), ])
p1counts <- dat[(nrow(dat) - 12 + 1):(nrow(dat) - 12 + 6), 1]
p1size <- rowSums(dat[(nrow(dat) - 12 + 1):(nrow(dat) - 12 + 6), ])
p2counts <- dat[(nrow(dat) - 12 + 7):nrow(dat), 1]
p2size <- rowSums(dat[(nrow(dat) - 12 + 7):nrow(dat), ])
ploidy <- 6
prior_vec <- rep(1/(ploidy + 1), length = ploidy + 1)
r1vec <- bin_post(ncounts = p1counts, ssize = p1size, prior = prior_vec)
r2vec <- bin_post(ncounts = p2counts, ssize = p2size, prior = prior_vec)
pk <- seq(0, ploidy) / ploidy ## the possible probabilities
pk <- (1 - seq_error) * pk + seq_error * (1 - pk)
pival <- 0.9
p1geno <- which.max(r1vec)
p2geno <- which.max(r2vec)
alpha <- 0.01
beta <- 0.01
dbinommat <- mapply(FUN = stats::dbinom, ocounts, osize, MoreArgs = list(prob = pk))
ldbinommat <- mapply(FUN = stats::dbinom, ocounts, osize, MoreArgs = list(prob = pk, log = TRUE))
qarray <- get_q_array(ploidy)
objout <- up_obj(pival = pival, p1geno = p1geno, p2geno = p2geno, alpha = alpha,
beta = beta, ocounts = ocounts, osize = osize,
dbinommat = dbinommat, qarray = qarray, r1vec = r1vec,
r2vec = r2vec)
itermax <- 100
objvec <- rep(NA, length = itermax)
for(iter_index in 1:itermax) {
objout_old <- objout
trash <- capture_output(fout <- up_fix(pival = pival, p1geno = p1geno, p2geno = p2geno, alpha = alpha,
beta = beta, ocounts = ocounts, osize = osize,
ldbinommat = ldbinommat, qarray = qarray, r1vec = r1vec,
r2vec = r2vec, update_beta = TRUE, update_pi = TRUE, update_geno = TRUE))
p1geno <- fout$p1geno
p2geno <- fout$p2geno
pival <- fout$pival
beta <- fout$beta
alpha <- fout$alpha
objout <- up_obj(pival = pival, p1geno = p1geno, p2geno = p2geno, alpha = alpha,
beta = beta, ocounts = ocounts, osize = osize,
dbinommat = dbinommat, qarray = qarray, r1vec = r1vec,
r2vec = r2vec)
objvec[iter_index] <- objout
expect_true(objout - objout_old > - 10 ^ -12)
}
r1vec <- rep(1 / (ploidy + 1), length = ploidy + 1)
r2vec <- rep(1 / (ploidy + 1), length = ploidy + 1)
trash <- capture_output(maxout <- updog_maximize(ocounts = ocounts, osize = osize, qarray = qarray,
r1vec = r1vec, r2vec = r2vec, pk = pk, pival = 0.99,
alpha = 0.1, beta = 0.1,
est_fudge = TRUE,
tol = 10 ^ -4, itermax = 1000,
update_geno = TRUE, update_pi = TRUE, update_beta = TRUE))
plot_geno(ocounts = ocounts, osize = osize, ploidy = ploidy, p1counts = p1counts,
p2counts = p2counts, p1size = p1size, p2size = p2size, ogeno = maxout$ogeno, prob_ok = maxout$theta)
plot(sort(apply(maxout$opostprob, 2, max)))
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.