Last updated: 2017-01-15

Code version: 94a3347615361216b39b8a9333bfa8354794e0ff

First, we load the necessary libraries.

library(REBayes)
## Loading required package: Matrix
library(ashr)

The NPMLE

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

Session information

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