orient-bias-filter: Orientation Bias Filter

# Explore the limitation of numeric integration required for the
# variant identification task using the empirical Bayes model.

library(devtools);
load_all("..");

# Set up a problem where the profile likelihood
# of phi diverges at 0

# This causes the standard integration function to crash.

n <- 200;
alpha_phi <- 0.1;
beta_phi <- 10;
phi <-  alpha_phi / (alpha_phi + beta_phi);
e.mean <- 20;
e.sd <- 2;
theta <- 0;

data <- simulate_orient_bias_read_stats(n, theta, phi, e.mean, e.sd)


eps <- 1e-6;

phis <- seq(eps, 1-eps, 0.001);
#phis <- exp(seq(-10-eps, -eps, 0.5));

# log likelihood with known theta
plot(phis, unlist(lapply(phis, function(phi) lp_bases_given_params(data, theta, phi))), type="l");

# likeihood with known theta
plot(phis, unlist(lapply(phis, function(phi) exp(lp_bases_given_params(data, theta, phi)))), type="l");

xs <- unlist(lapply(phis, function(phi) lp_bases_given_params(data, theta, phi)));
sum(exp(xs))
exp(log_sum_exp(xs))

# log density of phi
plot(phis, unlist(lapply(phis, function(phi) dbeta(phi, alpha_phi, beta_phi, log=TRUE))), type="l");

# density of phi
plot(phis, unlist(lapply(phis, function(phi) dbeta(phi, alpha_phi, beta_phi))), type="l");

# profile likelihood of phi with known theta
ys <- unlist(lapply(phis,
	function(phi) {
		exp(
			lp_bases_given_params(data, theta, phi) + 
			dbeta(phi, alpha_phi, beta_phi, log=TRUE) + 
			dbeta(theta, 1, 1, log=TRUE)
		)
	}));
plot(phis, ys, type = "l");
# expectation of phi under its profile likelihood
sum(ys * phis) / sum(ys)
abline(v=phi, col="red");


# maximum likelihood estimator of phi
# with naive implementation
phi_hat_mle_naive <- phis[which.max(unlist(lapply(phis, function(phi) lp_bases_given_params(data, theta, phi))))];
phi_hat_mle_naive

# MLE of phi
phi_hat_mle <- estimate_phi(data, theta, 0.5);
phi_hat_mle

# MLE of phi given inflated theta
estimate_phi(data, theta + 0.1, 0.5);


# marginal evidence for null model with theta = 0
# using standard integration
log(integrate_over_unit(
	function(phi) {
		exp(
			lp_bases_given_params(data, 0, phi) +
				dbeta(phi, alpha_phi, beta_phi, log=TRUE)
				#dbeta(theta, 1, 1, log=TRUE)
		)
	}
)$value)



thetas <- seq(eps, 1-eps, 0.01);
#thetas <- exp(seq(-10-eps, -eps, 0.5))

# profile likelihood of theta with phi = phi_hat
plot(thetas, unlist(lapply(thetas, function(theta) exp(lp_bases_given_params(data, theta, phi_hat)))), type="l");
abline(v=theta, col="red")

# profile likelihood of theta with phi = 0
plot(thetas, unlist(lapply(thetas, function(theta) exp(lp_bases_given_params(data, theta, 0)))), type="l");
abline(v=theta, col="red")

plot(thetas, unlist(lapply(thetas, function(theta) exp(lp_bases_given_params(data, theta, 0)))), type="l");
abline(v=theta, col="red")

lp_f <- function(phi, theta) {
	lp_bases_given_params(data, theta, phi) +
		dbeta(phi, alpha_phi, beta_phi, log=TRUE) +
		dbeta(theta, 1, 1, log=TRUE)
}

p_f <- function(phi, theta) {
	exp(
		lp_bases_given_params(data, theta, phi) +
			dbeta(phi, alpha_phi, beta_phi, log=TRUE) +
			dbeta(theta, 1, 1, log=TRUE)
	)
}

# problem with standard integration:
# when f rises quickly at a boundary, the
# integration function crashes
log(standard_integration(function(x) p_f(x, theta=0), lower=0, upper=1)$value)

# y is extremely skewed, with virtually all
# of its mass near 0.0
plot(phis, unlist(lapply(phis, p_f, theta=0)), type="l");
# log y space appears smoother, but there is
# still a discontinguity near 0.0
plot(phis, unlist(lapply(phis, lp_f, theta=0)), type="l");

# the simpler integration methods do not
# crash only because they do not go close
# enough to the boundary
# therefore, the results are sensitive to
# how close the grid points are placed near
# the boundaries

log_midpoint_integration(phis, lp_f, theta=0)
log_trapezoid_integration(phis, lp_f, theta=0)

log(midpoint_integration(phis, p_f, theta=0))

# Integrate in real space (instead of the constrained [0, 1] space) after
# logit transform would allow us to approach the boundary as close as possible
# However, at one point, numeric overflow occurs and NaN results

# take logit(phi) as input
lphis <- seq(-5, 5, 0.2);
log(midpoint_integration_logit_x(lphis, p_f, theta=0));

# take logit(phi) as input
lphis <- seq(-10, 10, 0.2);
log(midpoint_integration_logit_x(lphis, p_f, theta=0));

lphis <- seq(-20, 20, 0.2);
log(midpoint_integration_logit_x(lphis, p_f, theta=0));

lphis <- seq(-30, 30, 1);
log(midpoint_integration_logit_x(lphis, p_f, theta=0));

lphis <- seq(-40, 40, 1);
log(midpoint_integration_logit_x(lphis, p_f, theta=0));

cphis <- centered_grid(0);
log(midpoint_integration(cphis, p_f, theta=0))

abline(v=cphis, col="blue")

p_f(0, theta=0)
p_f(1e-10, theta=0)
p_f(1e-12, theta=0)
p_f(1e-6, theta=0)
p_f(1e-4, theta=0)

####

# Try models implemeneted with different numeric integration methods

f <- orient_bias_identify_bayes_model_standard(data, alpha_phi, beta_phi);
f
exp(diff(f))

f.m <- orient_bias_identify_bayes_model_midpoint(data, alpha_phi, beta_phi);
f.m
exp(diff(f.m))

# same results as orient_bias_identify_bayes_model_midpoint
#f.lm <- orient_bias_identify_bayes_model_log_midpoint(data, alpha_phi, beta_phi);
#f.lm
#exp(diff(f.lm))

f.mc <- orient_bias_identify_bayes_model_midpoint_cgrid(data, alpha_phi, beta_phi);
f.mc
exp(diff(f.mc))

f.ml <- orient_bias_identify_bayes_model_midpoint_logit_x(data, alpha_phi, beta_phi);
f.ml
exp(diff(f.ml))

# problematic approximation
f.l <- orient_bias_identify_bayes_model_large_approx(data, alpha_phi, beta_phi);
f.l
exp(diff(f.l))