Last updated: 2017-01-15
Code version: 94a3347615361216b39b8a9333bfa8354794e0ff
First, we load the necessary libraries.
library(REBayes)
## Loading required package: Matrix
library(ashr)
Here we show how we can use ash to (approximately) compute the unconstrained NPMLE - that is, estimate the underlying distribution g by maximizing the likelihood without the unimodal constraint. See Koenker and Mizera, JASA 2014, for background.
Within ash one can approximate the npmle as a mixture of uniforms on a dense grid of non-overlapping values - this results in a piecewise constant density with changes in the density only at the grid points.
The following example comes from the REBayes vignette by Koenker and Gu. The underlying g is a mixture of a point mass (weight 0.8) at 0 and a point mass (weight 0.2) at 2.
set.seed(102)
y <- c(rep(0,800),rnorm(200,2)) + rnorm(1000)
z <- GLmix(y)
Now we fit the NPMLE using ‘ashr’ and compare with the REBayes solution.
grid = seq(from = min(z$x),to = max(z$x),length = 1000)
k = length(grid)
y.ash.npmle =
ash(y,1,g = unimix(pi = rep(1/(k-1),(k-1)),a = grid[-k],b = grid[-1]),
method = "shrink")
# Plot the model fitted using GLmix against the fitted ash model.
plot(z$x,cumsum(z$y)/sum(z$y),col = 2,main = "Estimated cdf",type = "l")
lines(cdf.ash(y.ash.npmle,x = z$x),type = "l",lty = "dashed",lwd = 2)
# Compare the likelihood given the GLmix model and given ash model.
print(z$logLik,digits = 12)
[1] -1643.86778293
print(y.ash.npmle$loglik,digits = 12)
[1] -1643.85995729
sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.2
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] knitr_1.15.1 ashr_2.0.5 REBayes_0.63 Matrix_1.2-7.1
loaded via a namespace (and not attached):
[1] Rcpp_0.12.8 magrittr_1.5 MASS_7.3-45
[4] doParallel_1.0.10 pscl_1.4.9 SQUAREM_2016.8-2
[7] lattice_0.20-34 foreach_1.4.3 stringr_1.1.0
[10] tools_3.3.2 parallel_3.3.2 grid_3.3.2
[13] htmltools_0.3.5 iterators_1.0.8 assertthat_0.1
[16] yaml_2.1.14 rprojroot_1.1 digest_0.6.10
[19] codetools_0.2-15 evaluate_0.10 rmarkdown_1.3
[22] stringi_1.1.2 Rmosek_7.1.3 backports_1.0.4
[25] truncnorm_1.0-7