Compare_ebpmf_nmf

Goal

In previous experiments https://zihao12.github.io/ebpmf_demo/Issue_ebpmf_issue2.html, we found ebpmf and nmf (I use Lee’s multiplicative update, so it is EM algorithm for solving the MLE) have different performance on different dataset. Here I want to simulate data from a quantifiable noise level, as described here https://www.overleaf.com/project/5bd084d90a33772e7a7f99a2 (pve.tex), and see their performance.

Data setup and Comparison Metric

\[ \begin{align} & X_{ij} \sim Pois(\Lambda_{ij} U_{ij})\\ & \Lambda_{ij} = \sum_k L_{ik} F_{jk}\\ & L_{ik} \sim Gamma(a_L, b_L)\\ & F_{jk} \sim Gamma(a_F, b_F)\\ & U_{ij} \sim Gamma(1/a, 1/a)\\ \end{align} \]

Given \(a_L, a_F, b_L, b_F\), I simulate \(L, F\), and get \(\Lambda_{true} = L F^t\). \ Then for a chosen level of pve (according to the writeup), I compute the corresponding \(a\), then sample \(Y_{train}\), and \(Y_{val}\) by: firstsampling \(U_{train}, U_{val}\), then forming the mean matrix and sampling from the Poisson.\ “ll_train” is \(log p(Y_{train}|\hat{\Lambda})\); “ll_val” is \(log p(Y_{val}|\hat{\Lambda})\); “RMSE” is \(RMSE(\Lambda_{true}, \hat{\Lambda})\).

simulate_data <- function(n, p, K, params, seed = 123){
  set.seed(seed)
  L = matrix(rgamma(n = n*K, shape = params$al, rate = params$bl), ncol = K)
  F = matrix(rgamma(n = p*K, shape = params$af, rate = params$bf), ncol = K)
  Lam =  L %*% t(F)
  X = matrix(rpois(n*p, Lam), nrow = n)
  Y = matrix(rpois(n*p, Lam), nrow = n)
  pve = compute_pve(K, a =  params$a, al = params$al, af = params$af,bl = params$bl, bf = params$bf)
  return(list(params = params,Lam = Lam,X = X, Y = Y,  pve =  pve))
}


simulate_Lam <- function(n, p, K, params, seed = 123){
  set.seed(seed)
  L = matrix(rgamma(n = n*K, shape = params$al, rate = params$bl), ncol = K)
  F = matrix(rgamma(n = p*K, shape = params$af, rate = params$bf), ncol = K)
  Lam =  L %*% t(F)
  return(list(params = params,Lam = Lam,L = L, F = F))
}

simulate_X <- function(sim, params){
  n = nrow(sim$Lam)
  p = ncol(sim$Lam)
  U1 = matrix(rgamma(n*p, shape = 1/params$a, rate = 1/params$a), nrow = n)
  U_Lam1 = sim$Lam * U1
  X = matrix(rpois(n*p, U_Lam1), nrow = n)
  U2 = matrix(rgamma(n*p, shape = 1/params$a, rate = 1/params$a), nrow = n)
  U_Lam2 = sim$Lam * U2
  Y = matrix(rpois(n*p, U_Lam2), nrow = n) ## do we need to sample a different U?
  pve = compute_pve(K, a =  params$a, al = params$al, af = params$af,bl = params$bl, bf = params$bf)
  return(list(params = params,Lam = sim$Lam,X = X, Y = Y,  pve =  pve))
}


compute_pve <- function(K, a, al, af, bl, bf){
  var_lam = K * (1 + al + af) * al * af/(bl * bf)^2
  var_x = (a + 1) * var_lam + a * (K * al * af / (bl *  bf))^2
  return(var_lam/var_x)
}

compute_rmse <- function(lam1, lam2){
  return(sqrt(mean((lam1 - lam2)^2)))
}

compute_ll <- function(X, lam){
  return(sum(dpois(X,lam, log = T)))
}

## Fix  other params and choose "a" that can achieve the  targeted  pve
adjust_a_by_pve  <- function(params, K, pve){
  al = params$al
  af = params$af
  bl = params$bl
  bf = params$bf
  var_lam = K * (1 + al + af) * al * af/(bl * bf)^2
  numer = var_lam/pve - var_lam
  denom = var_lam + (K* al*af/(bl*bf))^2
  params[["a"]] = numer/denom
  return(params)
}

n = 100
p = 200
K = 2
params = list(al = 10, bl =  10, af = 10, bf = 10, a = NA)
m = 2  ## for ebpmf_exp
maxiter = 100

pve_ = 0.1*seq(1, 10, 2)
pves = c(); methods = c(); ll_trains  = c(); ll_vals  = c(); RMSEs = c();
sim_ = simulate_Lam(n, p, K, params)

for(pve in pve_){
  params = adjust_a_by_pve(params, K, pve)
  sim = simulate_X(sim_, params)
  out_ebpmf_exp = ebpmf::ebpmf_exponential_mixture(sim$X, K = K, m = m, maxiter.out = maxiter)
  out_ebpmf_exp[["lam"]] = out_ebpmf_exp$qg$qls_mean %*% t(out_ebpmf_exp$qg$qfs_mean)
  ll_trains = c(ll_trains, compute_ll(sim$X, out_ebpmf_exp[["lam"]]))
  ll_vals = c(ll_vals, compute_ll(sim$Y,out_ebpmf_exp[["lam"]]))
  RMSEs = c(RMSEs, compute_rmse(sim$Lam, out_ebpmf_exp[["lam"]]))
  methods = c(methods,"ebpmf_exp")
  pves = c(pves, sim$pve)

  out_ebpmf_point = ebpmf::ebpmf_point_gamma(sim$X, K = K, maxiter.out = maxiter)
  out_ebpmf_point[["lam"]] = out_ebpmf_point$qg$qls_mean %*% t(out_ebpmf_point$qg$qfs_mean)
  ll_trains = c(ll_trains, compute_ll(sim$X, out_ebpmf_point[["lam"]]))
  ll_vals = c(ll_vals, compute_ll(sim$Y,out_ebpmf_point[["lam"]]))
  RMSEs = c(RMSEs, compute_rmse(sim$Lam, out_ebpmf_point[["lam"]]))
  methods = c(methods,"ebpmf_point")
  pves = c(pves, sim$pve)

  out_nnmf = NNLM::nnmf(sim$X, k = K, max.iter = maxiter, method = "lee")
  out_nnmf[["lam"]] = out_nnmf$W %*% out_nnmf$H
  ll_trains = c(ll_trains, compute_ll(sim$X, out_nnmf[["lam"]]))
  ll_vals = c(ll_vals, compute_ll(sim$Y,out_nnmf[["lam"]]))
  RMSEs = c(RMSEs, compute_rmse(sim$Lam, out_nnmf[["lam"]]))
  methods = c(methods,"nnmf")
  pves = c(pves, sim$pve)
}

df = data.frame(pve = pves, method = methods, ll_train = ll_trains, ll_val = ll_vals, RMSE = RMSEs)

Below are the results (the black line in the first two plots are likelihood for oracle)

ggplot(df)+
  geom_point(aes(x = pve, y = ll_train, color = method))+
  geom_line(aes(x = pve, y = ll_train, color = method))+
  geom_abline(slope = 0, intercept = compute_ll(sim$X, sim$Lam))

Past versions of plot-1.png

ggplot(df)+
  geom_point(aes(x = pve, y = ll_val, color = method))+
  geom_line(aes(x = pve, y = ll_val, color = method))+
  geom_abline(slope = 0, intercept = compute_ll(sim$Y, sim$Lam))

Past versions of plot-2.png

ggplot(df)+
  geom_point(aes(x = pve, y = RMSE, color = method)) +
  geom_line(aes(x = pve, y = RMSE, color = method))

Past versions of plot-3.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] ggplot2_3.2.1  NNLM_0.4.2     gtools_3.8.1   ebpmf_0.1.0   
[5] testthat_2.2.1

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