fix converting BF to posterior

NAMESPACE

@@ -11,6 +11,7 @@ S3method(fgem_lik,dgCMatrix)
 S3method(fgem_lik,matrix)
 export("%>%")
 export("%||%")
+export(BF2posterior)
 export(cv_fgem)
 export(cv_relax_fgem)
 export(fgem)

R/FGEM_Func.R

@@ -491,7 +491,7 @@ fgem_null_lik <- function(BF,log_BF=FALSE) {
 
 fgem_null_fit <- function(BF, log_BF = FALSE, ...) {
     pta <- proc.time()
-    lbf <- stats::optimize(fgem_lik_stan,
+    lbf <- stats::optimize(fgem_lik,
                     lower = -10,
                     upper = 100, X = matrix(0, nrow = length(BF), ncol = 0), BF = BF,
                     log_BF=log_BF,

R/utils.R

@@ -34,7 +34,8 @@ summarise_cv_lik <- function(fit, summarise_Beta = TRUE) {
                        cv_sum = sum(cv_lik),
                        mean_time=mean(time),
                        l0_mean = mean(l0n),
-                       l0_sd = sd(l0n),
+                       l1_mean = mean(l1n),
+                       l2_mean = mean(l2n)
                        ) %>%
             dplyr::group_by(!!!mgrps)
     }else{
@@ -44,8 +45,10 @@ summarise_cv_lik <- function(fit, summarise_Beta = TRUE) {
                        Beta = summ_beta(Beta,
                                         drop_0 = FALSE),
                        cv_sum = sum(cv_lik),
+                       mean_time=mean(time),
                        l0_mean = mean(l0n),
-                       l0_sd = sd(l0n)
+                       l1_mean = mean(l1n),
+                       l2_mean = mean(l2n)
                    ) %>%
             dplyr::group_by(!!!mgrps)
     }
@@ -318,22 +321,72 @@ empty_x <- function(x, rownames = names(x)) {
         matrix(0, nrow = nr, ncol = 0, dimnames = list(rownames, NULL))
 }
 
-predict_uniform_prior <- function(BF, Genes = names(BF)) {
-    predict_fgem(fgem_null_fit(BF), empty_x(BF, Genes))
+predict_uniform_prior <- function(BF, Genes = names(BF), log_BF = FALSE) {
+        predict_fgem(fgem_null_fit(BF, log_BF = log_BF), empty_x(BF, Genes), log = log_BF)
 }
 
-BF2posterior <- function(BF,Genes=names(BF), prior = predict_uniform_prior(BF,Genes)) {
-    return((prior * BF) / ((prior * BF) + (1 - prior)))
+##' @title convert BF to posterior estimates
+##'
+##'
+##' @param BF bayes factor (or log bayes factor
+##' @param Genes character vector for gene names
+##' @param prior prior (or log prior)
+##' @param log_BF whether to compute in log space (and return in log space)
+##' @return
+##' @export
+BF2posterior <- function(BF,Genes=names(BF), prior = predict_uniform_prior(BF,Genes),log_BF=FALSE) {
+    if(!log_BF)
+        return((prior * BF) / ((prior * BF) + (1 - prior)))
+    return(log_BF2_logpost(BF, stats::qlogis(prior, log.p = TRUE)))
 }
 
 
 log_BF2_logpost <- function(logBF, xb) {
-        apply(xb - logBF, 2, log_1p_exp)
+  -log1pexp(xb + logBF) + xb + logBF
+  ## #    apply(xb - logBF, 2, log_1p_exp)
 }
 
 
 
 
+##' @title Compute  f(a) = log(1 - exp(-a))  stably
+##' @param a numeric vector of positive values
+##' @param cutoff  log(2) is optimal, see  Maechler (201x) .....
+##' @return f(a) == log(1 - exp(-a)) == log1p(-exp(-a)) == log(-expm1(-a))
+##' @author Martin Maechler, May 2002 .. Aug. 2011
+log1mexp <- function(a, cutoff = log(2)) ## << log(2) is optimal >>
+{
+    if(has.na <- any(ina <- is.na(a))) {
+	y <- a
+	a <- a[ok <- !ina]
+    }
+    if(any(a < 0))## a == 0  -->  -Inf	(in both cases)
+	warning("'a' >= 0 needed")
+    tst <- a <= cutoff
+    r <- a
+    r[ tst] <- log(-expm1(-a[ tst]))
+    r[!tst] <- log1p(-exp(-a[!tst]))
+    if(has.na) { y[ok] <- r ; y } else r
+}
+
+##' @title Compute  f(x) = log(1 + exp(x))  stably and quickly
+log1pexp <- function(x, c0 = -37, c1 = 18, c2 = 33.3)
+{
+    if(has.na <- any(ina <- is.na(x))) {
+	y <- x
+	x <- x[ok <- !ina]
+    }
+    r <- exp(x)
+    if(any(i <- c0 < x & (i1 <- x <= c1)))
+	r[i] <- log1p(r[i])
+    if(any(i <- !i1 & (i2 <- x <= c2)))
+	r[i] <- x[i] + 1/r[i] # 1/exp(x) = exp(-x)
+    if(any(i3 <- !i2))
+	r[i3] <- x[i3]
+    if(has.na) { y[ok] <- r ; y } else r
+}
+
+
 ##' @title Convert a table-style sparse matrix representation to a sparse matrix
 ##'
 ##'

tests/testthat/testEM.R

@@ -19,84 +19,45 @@ test_that("FGEM works on simulated data",{
   bf <- log(pxz1)-log(pxz2)
   ebf <- exp(bf)
 
-#   new_grad <- function(par,X,eBF){
-#       myenv <- new.env()
-#       assign("Beta",par, envir = myenv)
-# #      assign("Beta0",par[1] , envir = myenv)
-#       assign("eBF", eBF, envir = myenv)
-#       assign("X", X, envir = myenv)
-#       rd <- numericDeriv(quote(sum(log(eBF + (1 - eBF)/t(1 + exp(X%*%Beta[-1]+Beta[1]))))),c("Beta"),myenv)
-#       return(c(attr(rd,"gradient")))
-#   }
-  plan(multisession)
-  tr <- fgem:::fgem_bfgs(X,bf,log_BF=TRUE ,alpha=.5,lambda=lambda,hess=TRUE,grad=TRUE,verbose=FALSE)
-  ctr <- fgem::cv_fgem(X,ebf,log_BF=FALSE ,alpha=.5,hess=TRUE,grad=TRUE,verbose=FALSE,nlambda=4,v=4)
-  rlct <- cv_relax_fgem(X,ebf,log_BF=FALSE ,alpha=.5,hess=TRUE,grad=TRUE,verbose=FALSE,nlambda=4,v=4)
-  sctr <- fgem:::summarise_cv_lik(ctr)
-  check_lik <- purrr::pmap_dbl(tr,function(Beta,l1,l2,...){
-    B <- Beta$Beta
-    l <- -fgem::fgem_lik(B,X,bf,log_BF=TRUE)
-    l1p <- fgem:::l1_norm_penalty(B,l1 = l1)
-    l2p <- fgem:::l2_norm_penalty(B,l2 = l2)
-    return(-(l+l1p+l2p))
-  })
+  rlct <- fgem::cv_relax_fgem(X,ebf,log_BF=FALSE ,alpha=.5,hess=TRUE,grad=TRUE,verbose=FALSE,nlambda=4,v=4)
+  sctr <- fgem:::summarise_cv_lik(rlct)
 
-  tB0 <- tr$Beta[[1]]$Beta
-  check_grad <- purrr::pmap(tr,function(l1,l2,...){
-    g <- fgem::fgem_grad(tB0,X,bf,log_BF=TRUE,l1=l1,l2=0)
-    return(g)
-  })
-  gmat0 <- do.call(cbind,check_grad)
+})
+
+
+testthat::test_that("uninformative BF + uninformative prior -> uninformative posterior",{
 
-  rgrd <- map(fgem::fgem_grad(par = tr$Beta[[1]]$Beta,X,bf,l2=0,l1=tr$l1[1],log_BF=TRUE)[-1])
+  BF <- 1
+  hBF <- BF + 0.1
+  lBF <- BF - 0.1
 
+  prior <- 0.5  
+  hprior <- prior + 0.1
+  lprior <- prior - 0.1
 
 
-  mutate(tr,check_lik=check_lik) %>% filter(convergence==0,lik<max(lik)) %>% ggplot(aes(x=lik,y=check_lik))+geom_point()
-  fgem:::fgem_lik(tr$Beta[[3]]$Beta,X,bf,0,0,log_BF=TRUE)
-  new_lik(tr$Beta[[3]]$Beta,X,exp(bf))
-  tg <- fgem:::fgem_grad(tr$Beta[[3]]$Beta,X,bf,0,0,log_BF=TRUE)
-   ebf <- exp(bf)
-  cg <- new_grad(par,X,eBF)
-
-  mb <- microbenchmark::microbenchmark(  nl = fgem:::fgem_lik(tr$Beta[[3]]$Beta,X,ebf,0,0,log_BF=FALSE),
-                                         ll = fgem:::fgem_lik(tr$Beta[[3]]$Beta,X,bf,0,0,log_BF=TRUE))
-
+  balanced_post <- fgem:::BF2posterior(BF,prior=prior,log_BF=FALSE)
+  expect_equal(balanced_post,0.5)
 
-  wl <- which(tr$l0n>1)[1]
-  l <- lambda[wl]
-  prec <- tr$l2[wl]
-  l1 <- tr$l1[wl]
-  l <- lambda[wl]
-  fgl <- fgem:::fgem_elasticnet(X,bf,alpha=0.5,lambda=l,log_BF=TRUE,grad=TRUE,verbose=TRUE)
+  hbhp <- fgem:::BF2posterior(hBF,prior=hprior,log_BF=FALSE)
+  expect_gt(hbhp,balanced_post)
 
+  lblp <- fgem:::BF2posterior(lBF,prior=lprior,log_BF=FALSE)
+  expect_lt(lblp,balanced_post)
 
-
-
-  gr <- fgem::fgem_grad(fgl$Beta[[1]]$Beta,X,BF=bf,log_BF=TRUE,prec=fgl$l2,l1=fgl$l1)
-  sctg <- sign(gr)
-  i <- 4
-  par <- fgl$Beta[[1]]$Beta
-  prec <- fgl$l2
-  l1 <- fgl$l1
-  BF <- bf
-  log_BF <- TRUE
-  lik_idx <- function(i,par,X,BF,log_BF,prec,l1,sdd=2){
-    par[i] <- rnorm(1,mean=par[i],sd=sdd)
-    nlik <- fgem_lik(par,X,BF,prec=prec,log_BF=log_BF,l1=l1)
-    uplik <-fgem_lik(par,X,BF,prec=0.0,log_BF=log_BF,l1=0.0)
-    gr <- fgem_grad(par,X,BF,prec=prec,log_BF=log_BF,l1=l1)
-    upgr <- fgem_grad(par,X,BF,prec=0.0,log_BF=log_BF,l1=0.0)
-    return(tibble::tibble(par=par[i],lik=nlik,uplik=uplik,grad=gr[i],upgrad=upgr[i]))
-  }
-  table(y[X[,1]==1])
-  tdf <- purrr::map_dfr(rep(2,1000),lik_idx,par=par,X=X,BF=BF,log_BF=log_BF,prec=prec,l1=l1,sdd=.5)
-  ggplot2::ggplot(tdf,ggplot2::aes(x=par,y=lik,col=grad))+ggplot2::geom_point()+ggplot2::geom_vline(xintercept=fgl$Beta[[1]]$Beta[2])+ggplot2::scale_color_gradient2()
-  ggplot(tdf,aes(x=par,y=uplik,col=upgrad))+geom_point()+geom_vline(xintercept=fgl$Beta[[1]]$Beta[4])+scale_color_gradient2()
 
-
+  lbalanced_post <- fgem:::BF2posterior(log(BF),prior=log(prior),log_BF=TRUE)
+  expect_equal(exp(lbalanced_post),0.5)
+
+
+  lhbhp <- exp(fgem:::BF2posterior(log(hBF),prior=log(hprior),log_BF=TRUE))
+  expect_equal(hbhp,lhbhp)
+
+
+  llblp <- exp(fgem:::BF2posterior(log(lBF),prior=log(lprior),log_BF=TRUE))
+  expect_equal(lblp,llblp)
+
 })
 
 
 
-