Issues_ebpmf_rankk

Experiment setup

• I put my code for https://zihao12.github.io/ebpmf_demo/ebpmf_rankk_demo2.html into a pacakge: https://github.com/zihao12/ebpmf

• It is the same as in https://zihao12.github.io/ebpmf_demo/ebpmf_rankk_demo2.html except that I record the ELBO function (from ebpmf::ebpmf_exponential_mixture_experiment function) and find a few issues.

Issues:

• ELBO does not strictly increase. I think, if we can really solve each ebpm problem optimally within that prior family \(\mathcal{G}\), we should mamke ELBO nondecreasing (which is why coordinate ascend works?). But since in practice, we only optimize it over a subset of \(\mathcal{G}\) (based on our chosen grid), we may not solve each subproblem well.  

• [solved] Numerical issue when computing L: at iteration 107, I get Error in verify.likelihood.matrix(L) : Input argument "L" should be a numeric matrix with >= 2 columns, >= 1 rows, all its entries should be non-negative, finite and not NA, and some entries should be positive. It happens in ebpm when computing \(l_{ik} = log \ NB(x_i; a_k, \frac{b_k}{b_k + s})\). Only the first column of l (corresponding to the smallest chosen mean (biggest rate)) are all -inf. Why it happens? This is how we choose the grid

  xprime = x
  xprime[x == 0] = xprime[x == 0] + 1
  mu_grid_min =  0.05*min(xprime/s)
  mu_grid_max = 2*max(x/s)

  In this case, the smallest nonzero \(x\) is very small, so mu_grid_min becomes very small (\(10^{-18}\)). As a result, the loglikelihood is -Inf for nonzero \(x\) (which is fine, as exp(-Inf) = 0), and NAN for zeros. The latter is strange as my code does not handle big - small number very well, and R produces NAN for 0*-Inf:

  x = 0; prob = 1 - 1e-20;
  x*log(1-prob)
  ## [1] NaN 

  I use this for my implementation of the “continous” negative binomial loglikelihood. I solve it (probably not the best way) by setting that term to 0 “manually”.

I have solved issue 2 and reproduce issue 1 here.
## Experiments

library(NNLM)
library(gtools)

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

## simulate a poisson mean problem
## to do:
simulate_pm  <-  function(n, p, dl, df, K,scale_b = 10, seed = 123){
  set.seed(seed)
  ## simulate L
  a = replicate(dl,1)
  b = 10*runif(dl)
  pi <- rdirichlet(1,rep(1/dl, dl))
  gl = list(pi = pi, a = a, b= b)
  L = matrix(replicate(n*K, sim_mgamma(gl)), ncol = K)
  ## simulate F
  a = replicate(df,1)
  b = 10*runif(df)
  pi <- rdirichlet(1,rep(1/df, df))
  gf = list(pi = pi, a = a, b= b)
  F = matrix(replicate(p*K, sim_mgamma(gf)), ncol = K)
  ## simulate X
  lam = L %*% t(F)
  X = matrix(rpois(n*p, lam), nrow = n)
  Y = matrix(rpois(n*p, lam), nrow = n)
  ## prepare output
  g = list(gl = gl, gf = gf)
  out = list(X = X, Y = Y, L = L, F = F, g = g)
  return(out)
}

Simulate data

### simulate data
n = 100
p = 200
K = 2
dl = 10
df = 10
scale_b = 5
sim = simulate_pm(n, p, dl, df, K, scale_b = scale_b)

Run ebpmf and plot ELBO.

## ebpmf
out = ebpmf::ebpmf_exponential_mixture_experiment(sim$X, K, maxiter.out = 150, verbose = F)

[1] "summary of  runtime:"
[1] "init           : 0.029000"
[1] "Ez     per time: 0.001277"
[1] "rank1  per time: 0.025443"

elbo_ebpmf= out$ELBO
out_ebpmf = out$qg
distance_max =  max(elbo_ebpmf) - elbo_ebpmf
plot(1:length(distance_max), distance_max, xlab = "niter", ylab = "distance to max (ELBO)", type = "l", col = "blue")

Past versions of unnamed-chunk-3-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] gtools_3.8.1 NNLM_0.4.2  

loaded via a namespace (and not attached):
 [1] workflowr_1.4.0 Rcpp_1.0.2      digest_0.6.21   rprojroot_1.3-2
 [5] ebpmf_0.1.0     backports_1.1.5 git2r_0.25.2    magrittr_1.5   
 [9] evaluate_0.14   ebpm_0.0.0.9000 stringi_1.4.3   fs_1.3.1       
[13] whisker_0.3-2   rmarkdown_1.13  tools_3.5.1     stringr_1.4.0  
[17] glue_1.3.1      mixsqp_0.1-120  xfun_0.8        yaml_2.2.0     
[21] compiler_3.5.1  htmltools_0.3.6 knitr_1.25