Using ash for NPMLE

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