Last updated: 2017-01-15

Code version: 086c3157516dcce3907fd3fd6857a0aec1f35bfa

Outline

Ash has essentially three stages: - compute the likelihood matrix (whose j,k element is p(betahat_j| component k)) - estimate g - compute posterior quantities like E(beta_j | betahat_j, ghat)

Here we start to assess how long each step takes

First for normalmix:

library(ashr)
set.seed(100)
z = rnorm(100000,0,2)

total_time= system.time(res <- ash(z,1,mixcompdist="normal"))
fitted_g  = get_fitted_g(res)
data = res$data

time_no_output= system.time(ash(z,1,mixcompdist="normal",outputlevel=0))
time_min_output= system.time(ash(z,1,mixcompdist="normal",outputlevel=1))
time_no_output_no_fit= system.time(ash(z,1,mixcompdist="normal",outputlevel=0,fixg=TRUE,g=fitted_g))

time_matrix = system.time(llik <- t(ashr:::log_comp_dens_conv.normalmix(fitted_g,data)))

total_time
   user  system elapsed 
 28.889   3.035  16.617 
time_matrix
   user  system elapsed 
  0.112   0.053   0.166 
time_no_output
   user  system elapsed 
 21.823   1.082   7.442 
time_min_output
   user  system elapsed 
 24.233   1.665   9.565 
time_no_output_no_fit
   user  system elapsed 
  0.003   0.000   0.003 

Conclusion: the initial matrix calculation is relatively negligible. The full output takes about half the elapsed time, and a quarter of the user time.

Same for unimix:

total_time= system.time(res <- ash(z,1,mixcompdist="uniform"))
fitted_g  = get_fitted_g(res)
data = res$data

time_no_output= system.time(ash(z,1,mixcompdist="uniform",outputlevel=0))
time_min_output= system.time(ash(z,1,mixcompdist="uniform",outputlevel=1))

time_no_output_no_fit= system.time(ash(z,1,mixcompdist="uniform",outputlevel=0,fixg=TRUE,g=fitted_g))
time_matrix = system.time(llik <- t(ashr:::log_comp_dens_conv.unimix(fitted_g,data)))

total_time
   user  system elapsed 
 38.354   5.229  30.234 
time_matrix
   user  system elapsed 
  0.393   0.091   0.491 
time_no_output
   user  system elapsed 
 17.827   1.056   6.739 
time_min_output
   user  system elapsed 
 27.379   3.309  17.706 
time_no_output_no_fit
   user  system elapsed 
  0.002   0.001   0.002 

Now the full output takes more than half the elapsed time, and about half the user time.

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  

This site was created with R Markdown