demo/simulate_identification.R

In djhshih/orient-bias-filter: Orientation Bias Filter

library(io);
library(devtools);

load_all("..")

out.fname <- filename("obrs", path="sim-data", date=NA);

set.seed(1);

btest_to_list <- function(b) {
	g <- b;
	g$method <- NULL;
	g
}

htest_to_list <- function(h) {
	g <- h;
	class(g) <- "list";
	g$method <- NULL;
	g$statistic <- as.list(g$statistic);
	g$parameter <- as.list(g$parameter);
	g$estimate <- as.list(g$estimate);
	g
}

params <- list(
	n = 200,
	e.mean = 30,
	e.sd = 2,
	theta = NA,
	phi = NA
);

# Estimate alpha parameter of Beta distribution based on moment matching.
estimate_alpha_match <- function(m, v2) {
	m * ( m * (1 - m) / v2 - 1)
}

# Estimate beta parameter of Beta distribution based on moment matching.
estimate_beta_match <- function(m, v2) {
	(1 - m) * ( m * (1 - m) / v2 - 1)
}

process <- function(params, attribs) {
	s <- with(params, simulate_orient_bias_read_stats(n, theta, phi, e.mean, e.sd));

	print(attribs)
	print(params)

	print(with(s, table(base0, base)))
	print(with(s, table(oriented, damage)))

	h.known.phi <- orient_bias_adjusted_variant_test_with_known_phi(s, params$phi);
	h.unknown.phi <- orient_bias_adjusted_variant_test(s, params$phi);

	# moment match with Var[phi] = sqrt( E[Phi] )
	alpha_phi <- estimate_alpha_match(params$phi, params$phi^2);
	beta_phi  <- estimate_beta_match(params$phi, params$phi^2);

	b.uncertain.phi <- orient_bias_adjusted_variant_test_with_variable_phi(s, alpha_phi = alpha_phi, beta_phi = beta_phi);

	print(h.known.phi)
	print(h.unknown.phi)
	print(b.uncertain.phi)

	qwrite(params, insert(out.fname, tag = c(attribs, "params"), ext = "yaml"));

	qwrite(s, insert(out.fname, tag = c(attribs, "data"), ext = "tsv"));

	qwrite(htest_to_list(h.known.phi), insert(out.fname, tag = c(attribs, "results", "known-phi"), ext = "yaml"));
	qwrite(htest_to_list(h.unknown.phi), insert(out.fname, tag = c(attribs, "results", "unknown-phi"), ext = "yaml"));
	qwrite(btest_to_list(b.uncertain.phi), insert(out.fname, tag = c(attribs, "results", "uncertain-phi"), ext="yaml"));
}

params$theta <- 0;
params$phi <- 0.01;
attribs <- c("no-signal", "low-damage");
process(params, attribs);

params$theta <- 0.01;
params$phi <- 0.01;
attribs <- c("low-signal", "low-damage");
process(params, attribs);

params$theta <- 0.05;
params$phi <- 0.01;
attribs <- c("med-signal", "low-damage");
process(params, attribs);

params$theta <- 0.1;
params$phi <- 0.01;
attribs <- c("high-signal", "low-damage");
process(params, attribs);

params$theta <- 0.0;
params$phi <- 0.1;
attribs <- c("no-signal", "high-damage");
process(params, attribs);

params$theta <- 0.01;
params$phi <- 0.1;
attribs <- c("low-signal", "high-damage");
process(params, attribs);

params$theta <- 0.05;
params$phi <- 0.1;
attribs <- c("med-signal", "high-damage");
process(params, attribs);

params$theta <- 0.1;
params$phi <- 0.1;
attribs <- c("high-signal", "high-damage");
process(params, attribs);