This is a simple little library that uses optim and mvtnorm to estimate
likelihoods conditional on rejection.
Note that currently it _does not work_ for dimension > 2. You could
in principle use the simulated methods, but they are slow and unstable.
Suppose you have estimated a beta.hat, and is covariance is
beta.hat.sigma. Given this covariance matrix, you have selected
this beta hat because it is outside of the bounds
(-beta.hat.lower, beta.hat.upper), where these quantities can be
vectors, each applying to the corresponding element of beta hat.
Note that these bounds take into account
both the p-value of the test and the covariance matrix of your estimate.
To estimate the maximum likelihood of your truncated observation
and plot the likelihood, you can follow this example.
# Code example:
library(ggplot2)
source("multidimension_rejection_bias_correction_lib.R")
results <- NumericEstimateTruncatedMean(x=beta.hat,
sigma=beta.hat.sigma,
x.lower=beta.hat.lower,
x.upper=beta.hat.upper,
verbose=T)
# The output is the output of optim, so
# results$par contains the maximum likelihood estimate.
print(results$message)
print(results$par)
# Plot the likelihood if you like (this is the 1d case).
grid.values <- seq(-1, beta.hat * 1.5, length.out=500)
lik.grid <- sapply(grid.values,
function(x) {
NumericConditionalLogLik(x=beta.hat, beta=x,
x.lower=beta.hat.lower,
x.upper=beta.hat.upper,
sigma=beta.hat.sigma)
})
ggplot() + geom_line(aes(x=grid.values, y=exp(lik.grid))) +
geom_vline(aes(xintercept=beta.hat), col="blue") +
geom_vline(aes(xintercept=beta.hat.upper), color="gray", lwd=2) +
geom_vline(aes(xintercept=results$par), color="red")