Compare_ebvaepm_ebpm

show VAE results

vae_out = py_load_object("data/poisson-normal.pkl", pickle = "pickle")
vae_out = data.frame(vae_out)

## This is what data looks like
ggplot(vae_out)+
  geom_histogram(aes(x = x), bins = 100)

Past versions of unnamed-chunk-2-1.png

ggplot(vae_out)+
  geom_point(aes(x = x, y = posterior_vae))

Past versions of unnamed-chunk-2-2.png

Compare with ebpm_exponential_mixture

library(ebpm)
fit_ebpm_exp = ebpm_exponential_mixture(as.vector(vae_out$x), s = 1, m = 2^0.25)
vae_out[["posterior_ebpm_exp"]] = fit_ebpm_exp$posterior$mean
ggplot(vae_out)+
  geom_point(aes(x = x, y = posterior_ebpm_exp))

Past versions of unnamed-chunk-3-1.png

## biggest weight
max(fit_ebpm_exp$fitted_g$pi)

[1] 1

## the mean  of the exponential component corresponding to the biggest weight
fit_ebpm_exp$fitted_g$scale[which.max(fit_ebpm_exp$fitted_g$pi)]

[1] 10.76347

Compare with ebpm_point_gamma

library(ebpm)
fit_ebpm_point = ebpm_point_gamma(as.vector(vae_out$x), s = 1)
vae_out[["posterior_ebpm_point"]] = fit_ebpm_point$posterior$mean
ggplot(vae_out)+
  geom_point(aes(x = x, y = posterior_ebpm_point))

Past versions of unnamed-chunk-5-1.png

## fitted_g from point_gamma
class(fit_ebpm_point$fitted_g) = "data.frame"
fit_ebpm_point$fitted_g

         pi0    shape    scale
1 0.03153703 4.853747 2.156189

Compare RMSE

## rmse(fit_vae, lam)
sqrt(mean((vae_out$posterior_vae - vae_out$lam)^2))

[1] 2.63505

## rmse(fit_ebpm_exp, lam)
sqrt(mean((fit_ebpm_exp$posterior$mean - vae_out$lam)^2))

[1] 2.962624

## rmse(fit_ebpm_point, lam)
sqrt(mean((fit_ebpm_point$posterior$mean - vae_out$lam)^2))

[1] 2.645054

## rmse(mle, lam)
sqrt(mean((vae_out$x- vae_out$lam)^2))

[1] 3.197445

look at \(\hat{g}\)

sample_point_gamma_one <- function(point_gamma_){
  point_gamma_ = fit_ebpm_point$fitted_g
  if(rbinom(1,1, point_gamma_$pi0) == 1){
    return(0)
  }else{
    return(rgamma(1,shape = point_gamma_$shape, scale = point_gamma_$scale))
  }
}

sample_point_gamma  <- function(n, point_gamma_, seed = 123){
  set.seed(seed)
  out = replicate(n, sample_point_gamma_one(point_gamma_))
  return(out)
}

## simulate a poisson mean problem from mixture of exponential
sample_expmix  <-  function(n,gammamix_, seed = 123){
  set.seed(seed)
  a = gammamix_$shape
  b = 1/gammamix_$scale
  pi = gammamix_$pi
  lam = replicate(n, sim_mgamma(a, b, pi))
  return(lam)
}


sim_mgamma <- function(a,b,pi){
  idx = which(rmultinom(1,1,pi) == 1)
  return(rgamma(1, shape = a[idx], rate =  b[idx]))
}

n = length(vae_out$lam)
#hist(sample_point_gamma(n, fit_ebpm_point$fitted_g), breaks = 100, freq = F)

truth_df = data.frame(samples = vae_out$lam, method = "truth")
vae_df = data.frame(samples = vae_out$prior_vae, method = "vae")
point_df = data.frame(samples = sample_point_gamma(n, fit_ebpm_point$fitted_g), method = "point_gamma")
exp_df   = data.frame(samples = sample_expmix(n, fit_ebpm_exp$fitted_g), method = "exponential_mixture")

samples_df = rbind(truth_df, vae_df, point_df, exp_df)

ggplot(samples_df, aes(samples, fill = method)) + geom_density(alpha = 0.2)+ggtitle("compare g_hat")

Past versions of unnamed-chunk-8-1.png

Session information

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] reticulate_1.12 ggplot2_3.2.1   ebpm_0.0.0.9001 testthat_2.2.1 

loaded via a namespace (and not attached):
 [1] gtools_3.8.1        tidyselect_0.2.5    xfun_0.8           
 [4] remotes_2.1.0       purrr_0.3.2         lattice_0.20-38    
 [7] colorspace_1.4-1    usethis_1.5.1       htmltools_0.3.6    
[10] yaml_2.2.0          rlang_0.4.0         pkgbuild_1.0.3     
[13] mixsqp_0.1-121      pillar_1.4.2        glue_1.3.1         
[16] withr_2.1.2         sessioninfo_1.1.1   stringr_1.4.0      
[19] munsell_0.5.0       gtable_0.3.0        workflowr_1.4.0    
[22] devtools_2.2.1.9000 memoise_1.1.0       evaluate_0.14      
[25] labeling_0.3        knitr_1.25          callr_3.2.0        
[28] ps_1.3.0            Rcpp_1.0.2          backports_1.1.5    
[31] scales_1.0.0        desc_1.2.0          pkgload_1.0.2      
[34] jsonlite_1.6        fs_1.3.1            digest_0.6.22      
[37] stringi_1.4.3       processx_3.3.1      dplyr_0.8.1        
[40] rprojroot_1.3-2     grid_3.5.1          cli_1.1.0          
[43] tools_3.5.1         magrittr_1.5        lazyeval_0.2.2     
[46] tibble_2.1.3        crayon_1.3.4        whisker_0.3-2      
[49] pkgconfig_2.0.3     ellipsis_0.3.0      Matrix_1.2-17      
[52] prettyunits_1.0.2   assertthat_0.2.1    rmarkdown_1.13     
[55] rstudioapi_0.10     R6_2.4.0            git2r_0.25.2       
[58] compiler_3.5.1