Test the fix point function and EM algorithm

Outline

I want to try incorporating the multi-resolution ideas from fast_ash.cpp into the fix point function in MixSquarem.cpp.

library(ashr)
set.seed(100)
nsamp=100000
z = rnorm(nsamp,0,2)
#now sort z so that they are in order
z = z[order(abs(z))]

res <- ash(z,1,mixcompdist="normal",outputlevel=4)
lik = res$fit_details$matrix_lik
fitted_g = get_fitted_g(res)

# set up the initial value of pi as uniform
pi = rep(1, ncomp(fitted_g))
control=list(K = 1, method=3, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,tol=1.e-07, maxiter=5000, trace=FALSE)

prior = rep(1,ncomp(fitted_g))

Rcpp::sourceCpp('MixSquarem.cpp')
fixptfn_orig(pi,lik,prior)

$fixedpointvector
 [1] 0.06212279 0.06222939 0.06233550 0.06254627 0.06296192 0.06376958
 [7] 0.06528978 0.06795816 0.07196982 0.07621205 0.07742674 0.07240996
[13] 0.06165978 0.04867429 0.03656324 0.02668879 0.01918194

$objfn
[1] 57641.2

fixptfn(pi,lik,prior)

$fixedpointvector
 [1] 0.06212279 0.06222939 0.06233550 0.06254627 0.06296192 0.06376958
 [7] 0.06528978 0.06795816 0.07196982 0.07621205 0.07742674 0.07240996
[13] 0.06165978 0.04867429 0.03656324 0.02668879 0.01918194

$objfn
[1] 57641.2

fixptfn(pi,lik,prior,tol=1e-8)

$fixedpointvector
 [1] 0.06212278 0.06222938 0.06233549 0.06254626 0.06296192 0.06376957
 [7] 0.06528978 0.06795817 0.07196983 0.07621206 0.07742675 0.07240997
[13] 0.06165978 0.04867430 0.03656325 0.02668879 0.01918194

$objfn
[1] 57641.21

system.time(fixptfn(pi,lik,prior))

   user  system elapsed 
  0.403   0.013   0.420 

system.time(fixptfn(pi,lik,prior,tol=1e-8))

   user  system elapsed 
  0.016   0.002   0.019 

# This is closer to the usual initial value we use in ash
normalize=function(x){x/sum(x)}
pi = rep(1/nsamp, ncomp(fitted_g))
pi[1]=1
pi = normalize(pi)

control$multiscale_tol=0
system.time(EM0<-cxxMixSquarem(lik,prior,pi,control))

   user  system elapsed 
162.368   5.755 169.487 

control$multiscale_tol=1e-5
system.time(EM5<-cxxMixSquarem(lik,prior,pi,control))

   user  system elapsed 
  1.573   1.897   3.513 

EM0

$pihat
 [1] 1.987464e-14 0.000000e+00 6.812328e-18 8.659740e-20 8.775203e-18
 [6] 1.929568e-18 2.907978e-15 5.612077e-11 8.617571e-06 8.176804e-01
[11] 1.823109e-01 0.000000e+00 3.395400e-15 3.656962e-09 0.000000e+00
[16] 0.000000e+00 0.000000e+00

$B
[1] 36854.42

$niter
[1] 139

$converged
[1] TRUE

EM5

$pihat
 [1] 5.965870e-15 0.000000e+00 8.157919e-18 6.654677e-18 4.576339e-18
 [6] 4.929390e-18 1.072553e-15 2.880648e-11 6.173629e-06 8.174338e-01
[11] 1.825600e-01 0.000000e+00 0.000000e+00 2.696856e-09 0.000000e+00
[16] 0.000000e+00 0.000000e+00

$B
[1] 36850.64

$niter
[1] 140

$converged
[1] TRUE

Session information

R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.5 (El Capitan)

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] ashr_2.1        workflowr_0.2.0 rmarkdown_1.3  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.8       rstudioapi_0.6    knitr_1.15.1     
 [4] magrittr_1.5      REBayes_0.73      MASS_7.3-45      
 [7] doParallel_1.0.10 pscl_1.4.9        SQUAREM_2016.8-2 
[10] lattice_0.20-34   foreach_1.4.3     stringr_1.1.0    
[13] tools_3.3.1       parallel_3.3.1    grid_3.3.1       
[16] git2r_0.18.0      htmltools_0.3.5   iterators_1.0.8  
[19] assertthat_0.1    yaml_2.1.14       rprojroot_1.1    
[22] digest_0.6.10     Matrix_1.2-7.1    codetools_0.2-15 
[25] evaluate_0.10     stringi_1.1.2     Rmosek_7.1.2     
[28] backports_1.0.4   truncnorm_1.0-7