compare ebpm with gauss_gh

Introduction

I want to compare ebpm with the algorithm Gauss-HG propsed in paper Bayesian inference on quasi-sparse count data. Below I first copy from their analysis http://dattahub.github.io/GHstancodes , then compare ebpm with theirs.

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(ggplot2)
theme_set(theme_bw())
library(plyr)
library(dplyr)
library(reshape2)

Gauss-HG algorithm

# setup Stan Gauss-HG sampler
{
  library(plyr)
  library(rstan)
  library(parallel)
  library(rbenchmark)
  
  #set_cppo("fast")
  stan.gh.code = "
  data{
  int<lower=0> J;
  int<lower=0> Y[J];
  real<lower=0> alpha;
  real<lower=0> a;
  real<lower=0> b;
  real<lower=0> gamma;
  real<lower=0> phi;
  }
  parameters{
  real<lower=0,upper=1> kappa[J];
  real<lower=0> theta[J];
  }
  model{
  for(i in 1:J) {
  increment_log_prob((a-1)*log(kappa[i])+(b-1)*log(1-kappa[i])-gamma*log(1-phi*kappa[i]));
  theta[i] ~ gamma(a, kappa[i]/(1-kappa[i]));
  Y[i] ~ poisson(theta[i]);
  }
  }
  "
  stan.gh.fit = stan_model(model_code=stan.gh.code, model_name="GH")
}

DIAGNOSTIC(S) FROM PARSER:
Info: increment_log_prob(...); is deprecated and will be removed in the future.
  Use target += ...; instead.

simulation

stan.iters = 10000
n.chains = 2
seed.val = 786
set.seed(seed.val)

n = 200; w = 0.9
y = rep(0,n); idx = rep(1,n)
lambdasparse = rep(0,n)
for (i in 1:n)
{
  if(i<=round(n*w)){
    lambdasparse[i]<-0.1
    idx[i] <- 0}
  else {lambdasparse[i] <-10}}

y = rpois(n,lambdasparse); 
gamma = mean(kmeans(y,centers=2)$centers)
alpha = 0.01
a = 0.5; b = 0.5
gh.data = list('J'=n,'Y'=y, 'alpha' = alpha,'a' = a, 'b' = b, 'gamma' = gamma, 'phi' = 0.99)

fit with Gauss-HG

{
  gh.res = sampling(stan.gh.fit, 
                     data = gh.data, 
                     iter = stan.iters,
                     warmup = floor(stan.iters/2),
                     thin = 2,
                     pars = c('kappa','theta'),
                     init = 0,
                     seed = seed.val, 
                     chains = 1)
  
  gh.theta.smpls = extract(gh.res, pars=c('theta'), permuted=TRUE)[[1]]
  gh.kappa.smpls = extract(gh.res, pars=c('kappa'), permuted=TRUE)[[1]]
  gh.theta.mean = apply(gh.theta.smpls,2,mean)
  gh.kappa.mean = apply(gh.kappa.smpls,2,mean)
  
  gh.sample.data = melt(extract(gh.res, permuted=TRUE))
  colnames(gh.sample.data) = c("iteration", "component", "value", "variable")
  
  gh.sample.data= gh.sample.data %>%
    filter(variable %in% c("theta","kappa")) 

  gh.sample.data.2 = gh.sample.data %>% group_by(component, variable) %>%
  summarise(upper = quantile(value, prob=0.975), 
            lower = quantile(value, prob=0.225),
            middle = mean(value))
}

SAMPLING FOR MODEL 'GH' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00019 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.9 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 20.3628 seconds (Warm-up)
Chain 1:                21.8004 seconds (Sampling)
Chain 1:                42.1632 seconds (Total)
Chain 1: 

fit with ebpm

library(ebpm)
fit_ebpm_gammamix = ebpm_gamma_mixture_single_scale(x = y, s = 1)
fit_ebpm_expmix = ebpm_exponential_mixture(x = y, s = 1)
fit_ebpm_pg = ebpm_point_gamma(x = y, s = 1)
fit_ebpm_tg = ebpm_two_gamma(x = y, s = 1)

fit_df = data.frame(
  data = y,
  lam_true = lambdasparse,
  gh = gh.theta.mean,
  ebpm_pg = fit_ebpm_pg$posterior$mean,
  ebpm_tg = fit_ebpm_tg$posterior$mean,
  ebpm_expmix = fit_ebpm_expmix$posterior$mean,
  ebpm_gammamix = fit_ebpm_gammamix$posterior$mean
)

ggplot(data = fit_df)+
  geom_point(aes(x = data, y = lam_true, color = "lam_true"))+
  geom_point(aes(x = data, y = gh, color = "gauss-hg"))+
  geom_point(aes(x = data, y = ebpm_pg, color = "ebpm_pg"))+
  geom_point(aes(x = data, y = ebpm_tg, color = "ebpm_tg"))

Take a closer look at those quasi-zeros (counts that comes from small lambda)

fit_df_small = fit_df[fit_df$lam_true < 1, ]
ggplot(data = fit_df_small)+
  geom_point(aes(x = data, y = lam_true, color = "lam_true"))+
  geom_point(aes(x = data, y = gh, color = "gauss_hg"))+
  geom_point(aes(x = data, y = ebpm_pg, color = "ebpm_point_gamma"))+
  geom_point(aes(x = data, y = ebpm_tg, color = "ebpm_two_gamma"))

Below I show the divergence between estimation and truth (Root Mean Squared Error, Kullback–Leibler divergence , Jensen-Shannon)

rmse <- function(true, est){
  return(sqrt(mean((true - est)^2)))
}
KL <- function(true,est){
  sum(ifelse(true==0,0,true * log(true/est)) + est - true)
}
JS  <- function(true,est){
  0.5*(KL(true, est) + KL(est, true))
}
RMSEs = c(rmse(lambdasparse, gh.theta.mean), rmse(lambdasparse, fit_ebpm_gammamix$posterior$mean), 
                     rmse(lambdasparse, fit_ebpm_expmix$posterior$mean), 
                     rmse(lambdasparse, fit_ebpm_pg$posterior$mean),
                     rmse(lambdasparse,fit_ebpm_tg$posterior$mean))

KLs = c(KL(lambdasparse, gh.theta.mean), KL(lambdasparse, fit_ebpm_gammamix$posterior$mean), 
                     KL(lambdasparse, fit_ebpm_expmix$posterior$mean), 
                     KL(lambdasparse, fit_ebpm_pg$posterior$mean),
                     KL(lambdasparse,fit_ebpm_tg$posterior$mean))

JSs = c(JS(lambdasparse, gh.theta.mean), rmse(lambdasparse, fit_ebpm_gammamix$posterior$mean), 
                     JS(lambdasparse, fit_ebpm_expmix$posterior$mean), 
                     JS(lambdasparse, fit_ebpm_pg$posterior$mean),
                     JS(lambdasparse,fit_ebpm_tg$posterior$mean))
data.frame(RMSE = RMSEs, KL = KLs, JS = JSs, row.names = c("guass-hg", "ebpm_gammamix", "ebpm_expmix", "ebpm_point_gamma", "ebpm_two_gamma"))

                      RMSE         KL         JS
guass-hg         1.1529090 88.3375455 60.4323648
ebpm_gammamix    0.9544262 24.7393189  0.9544262
ebpm_expmix      0.8374462 12.2507980 13.6031322
ebpm_point_gamma 0.9419141 24.2177962 32.0882989
ebpm_two_gamma   0.1140922  0.5261136  0.5411592

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

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] ebpm_0.0.0.9004    reshape2_1.4.3     dplyr_0.8.1       
[4] plyr_1.8.4         rstan_2.19.2       ggplot2_3.2.1     
[7] StanHeaders_2.19.0

loaded via a namespace (and not attached):
 [1] gtools_3.8.1       tidyselect_0.2.5   xfun_0.8          
 [4] purrr_0.3.2        colorspace_1.4-1   htmltools_0.3.6   
 [7] stats4_3.5.1       loo_2.1.0          yaml_2.2.0        
[10] rlang_0.4.1        pkgbuild_1.0.3     mixsqp_0.2-3      
[13] later_0.8.0        pillar_1.4.2       glue_1.3.1        
[16] withr_2.1.2        matrixStats_0.54.0 stringr_1.4.0     
[19] munsell_0.5.0      gtable_0.3.0       workflowr_1.5.0   
[22] evaluate_0.14      labeling_0.3       inline_0.3.15     
[25] knitr_1.25         callr_3.2.0        httpuv_1.5.1      
[28] ps_1.3.0           parallel_3.5.1     Rcpp_1.0.2        
[31] promises_1.0.1     scales_1.0.0       backports_1.1.5   
[34] fs_1.3.1           gridExtra_2.3      digest_0.6.22     
[37] stringi_1.4.3      processx_3.3.1     grid_3.5.1        
[40] rprojroot_1.3-2    cli_1.1.0          tools_3.5.1       
[43] magrittr_1.5       lazyeval_0.2.2     tibble_2.1.3      
[46] crayon_1.3.4       whisker_0.3-2      pkgconfig_2.0.3   
[49] prettyunits_1.0.2  assertthat_0.2.1   rmarkdown_1.13    
[52] R6_2.4.0           git2r_0.26.1       compiler_3.5.1