a

Author: william-denault

DESCRIPTION

@@ -23,7 +23,8 @@ Imports:
     ggplot2,
     cowplot(>= 1.1.3),
     smashr,
-    Rcpp
+    Rcpp,
+    ebnm,
 Suggests:
     susieR,
     knitr,

NAMESPACE

@@ -82,7 +82,6 @@ export(L_mixsq)
 export(L_mixsq.mixture_normal)
 export(L_mixsq.mixture_normal_per_scale)
 export(Pois_fSuSiE)
-export(Pois_fSuSiE2)
 export(TI_regression)
 export(TI_regression.susiF)
 export(affected_reg)
@@ -152,6 +151,7 @@ export(plot_susiF)
 export(plot_susiF_effect)
 export(plot_susiF_pip)
 export(pois_mean_GP)
+export(pois_smooth_split)
 export(post_mat_mean)
 export(post_mat_mean.mixture_normal)
 export(post_mat_mean.mixture_normal_per_scale)
@@ -161,6 +161,7 @@ export(post_mat_sd.mixture_normal_per_scale)
 export(remap_data)
 export(simu_IBSS_ash_vanilla)
 export(simu_IBSS_per_level)
+export(smash_dwt)
 export(smash_regression)
 export(smash_regression.susiF)
 export(susiF)
@@ -204,6 +205,11 @@ export(vga_pois_solver_newton_1iter)
 export(which_dummy_cs)
 export(which_dummy_cs.susiF)
 export(which_lowcount)
+<<<<<<< Updated upstream
+=======
+import(ebnm)
+import(wavethresh)
+>>>>>>> Stashed changes
 importFrom(ashr,ash)
 importFrom(ashr,calc_loglik)
 importFrom(ashr,cdf_post)

R/EB_poisson_mean_routines.R

@@ -469,3 +469,231 @@ pois_mean_GP_opt_obj_gradient = function(theta,x,s,beta,sigma2,n){
 pois_mean_GP_obj = function(x,s,beta,sigma2,m,v){
   return(sum(x*m-s*exp(m+v/2)-log(sigma2)/2-(m^2+v-2*m*beta+beta^2)/2/sigma2+log(v)/2))
 }
+
+
+
+
+
+
+
+#'@title Smooth over-dispersed Poisson sequence
+#'@param x data vector
+#'@param maxiter,tol max iteration and tolerance for stopping it.
+#'@param Eb_init,sigma2_init initial values of smooth mean and nugget effect.
+#'@examples
+#' set.seed(12345)
+#' n=2^9
+#' sigma=0.5
+#' mu=c(rep(0.3,n/4), rep(3, n/4), rep(10, n/4), rep(0.3, n/4))
+#' x = rpois(n,exp(log(mu)+rnorm(n,sd=sigma)))
+#' fit = pois_smooth_split(x,maxiter=30)
+#' plot(x,col='grey80')
+#' lines(exp(fit$Eb))
+#' fit$sigma2
+#' plot(fit$obj)
+#'@details The problem is
+#'\deqn{x_i\sim Poisson(\exp(\mu_i)),}
+#'\deqn{\mu_i\sim N(b_i,\sigma^2),}
+#'\deqn{\b_i\sim g(.).}
+#'@export
+
+pois_smooth_split = function(x,
+                             s = NULL,
+                             Eb_init = NULL,
+                             sigma2_init = NULL,
+                             est_sigma2 = TRUE,
+                             maxiter = 100,
+                             tol=1e-5,
+                             filter.number = 1,
+                             family = 'DaubExPhase',
+                             verbose=FALSE,
+                             printevery = 10,
+                             ebnm_params=list(mode=0),
+                             optim_method='L-BFGS-B'){
+
+  n = length(x)
+  if(is.null(s)){
+    s = 1
+  }
+  if(length(s)==1){
+    s = rep(s,n)
+  }
+  if(is.null(Eb_init)){
+    Eb = log(runmed(x/s,1 + 2 * min((n-1)%/% 2, ceiling(0.1*n)))+0.01)
+  }else{
+    Eb = Eb_init
+  }
+  if(is.null(sigma2_init)){
+    sigma2 = var(log(x/s+0.01)-Eb)
+  }else{
+    sigma2 = sigma2_init
+  }
+
+  W = (t(wavethresh::GenW(n,filter.number,family)))[-1,]
+
+  mu_pm = rep(0,n)
+  mu_pv = rep(1/n,n)
+  obj = -Inf
+
+  for(iter in 1:maxiter){
+    # get m, s^2
+    opt = optim(c(mu_pm,log(mu_pv)),
+                fn = pois_mean_GP_opt_obj,
+                gr = pois_mean_GP_opt_obj_gradient,
+                x=x,
+                s=s,
+                beta=Eb,
+                sigma2=sigma2,
+                n=n,
+                method = optim_method)
+    mu_pm = opt$par[1:n]
+    mu_pv = exp(opt$par[(n+1):(2*n)])
+    qb = smash_dwt(x=mu_pm,
+                   sigma=sqrt(sigma2),
+                   filter.number=filter.number,
+                   family=family,
+                   ebnm_params=ebnm_params,W=W)
+    Eb = qb$mu.est
+    Eb2 = qb$mu.est.var + Eb^2
+    # get sigma2
+    if(est_sigma2){
+      sigma2 = mean(mu_pm^2+mu_pv+Eb2-2*mu_pm*Eb)
+    }
+
+
+    # calc obj
+    obj[iter+1] = pois_smooth_split_obj(x,s,mu_pm,mu_pv,Eb,Eb2,sigma2,qb$dKL)
+
+    if(verbose){
+      if(iter%%printevery==0){
+        print(paste("Done iter",iter,"obj =",obj[iter+1]))
+      }
+
+    }
+
+    if((obj[iter+1]-obj[iter])<tol){
+      break
+    }
+
+  }
+  return(list(Emean = exp(mu_pm+mu_pv/2),
+              Vmean = exp(mu_pv-1)*exp(2*mu_pm+mu_pv),
+              Emu=mu_pm,
+              Vmu=mu_pv,
+              Eb=Eb,
+              Eb2=Eb2,
+              sigma2=sigma2,
+              obj=obj,
+              H = qb$dKL + sum(log(2*pi*mu_pv)/2-log(2*pi*sigma2)/2-(mu_pm^2+mu_pv-2*mu_pm*Eb+Eb2)/2/sigma2)))
+}
+
+pois_smooth_split_obj = function(x,s,m,s2,Eb,Eb2,sigma2,KLb){
+  return(sum(x*m-s*exp(m+s2/2)+log(s2)/2-log(sigma2)/2-(m^2+s2-2*m*Eb+Eb2)/2/sigma2)+KLb)
+}
+
+
+
+#'@title Empirical Bayes wavelet smoothing via DWT
+#'@description Smooth homogeneous Gaussian data.
+#'@param x data
+#'@param sigma known standard error
+#'@param filter.number,family wavelet family and filter number as in wavethresh package
+#'@param ebnm_params a list of `ebnm` parameters
+#'@param W the dwt matrix for calc posterior variance. Remove the first row which is all 1/sqrt(n)
+#'@return a list of
+#'  \item{mu.est:}{posterior mean}
+#'  \item{mu.est.var:}{posterior variance}
+#'  \item{loglik:}{log likelihood}
+#'  \item{dKL:}{KL divergence between g(the prior) and q(the posterior)}
+#'@import wavethresh
+#'@import ebnm
+#'@export
+smash_dwt = function(x,sigma,filter.number=1,
+                     family="DaubExPhase",
+                     ebnm_params=list(mode=0),W=NULL){
+
+  n = length(x)
+  J = log(n,2)
+  if(ceiling(J)!=floor(J)){
+    stop('Length of x must be power of 2')
+  }
+  if(is.null(ebnm_params)){
+    ebnm_params = ebnm_params_default()
+  }else{
+    temp = ebnm_params_default()
+    for(i in 1:length(ebnm_params)){
+      temp[[names(ebnm_params)[i]]] = ebnm_params[[i]]
+    }
+    ebnm_params = temp
+  }
+  tsum = sum(x)/sqrt(n)
+  x.w = wd(x, filter.number = filter.number,
+           family = family, type = "wavelet")
+
+  data.var = sigma^2
+  if(length(data.var==1)){
+    data.var = rep(data.var,n)
+  }
+
+  if(is.null(W)){
+    W = (t(GenW(n,filter.number,family)))[-1,]
+  }
+
+  if(length(sigma)==1){
+    x.w.v = rep(sigma^2,n-1)
+    tsum.var = sigma^2
+  }else{
+    x.w.v =  data.var
+    tsum.var = x.w.v[1]
+    x.w.v = x.w.v[-1]
+  }
+
+  dKL = 0
+  loglik.scale = c()
+  x.w.v.s = rep(0, 2^J-1)
+  for (j in 0:(J - 1)) {
+    x.pm = rep(0, 2^j)
+    #index = (((J - 1) - j) * n + 1):((J - j) * n)
+    index = (n-2^(j+1)+1):(n-2^j)
+    x.w.j = accessD(x.w, j)
+    x.w.v.j = x.w.v[index]
+    ind.nnull = (x.w.v.j != 0)
+
+    a = ebnm(x.w.j[ind.nnull],sqrt(x.w.v.j[ind.nnull]),
+             mode=ebnm_params$mode,
+             prior_family=ebnm_params$prior_family,
+             scale = ebnm_params$scale,
+             g_init = ebnm_params$g_init,
+             fix_g = ebnm_params$fix_g,
+             output = ebnm_params$output,
+             optmethod = ebnm_params$optmethod)
+
+    dKL = dKL + a$log_likelihood - Eloglik_GP(x.w.j[ind.nnull], sqrt(x.w.v.j[ind.nnull]),a$posterior$mean, a$posterior$mean^2+a$posterior$sd^2)
+    x.pm[ind.nnull] = a$posterior$mean
+    x.pm[!ind.nnull] = 0
+    x.w = putD(x.w, j, x.pm)
+    loglik.scale[j + 1] = a$log_likelihood
+    x.w.v.s[index[ind.nnull]] = a$posterior$sd^2
+    x.w.v.s[index[!ind.nnull]] = 0
+  }
+  mu.est = wr(x.w)
+  loglik = sum(loglik.scale)
+  #x.w.v.s = c(tsum.var,x.w.v.s)
+  mu.est.var = colSums(W^2*x.w.v.s)
+  return(list(mu.est=mu.est,mu.est.var=mu.est.var,loglik = loglik,dKL = dKL))
+}
+
+
+
+Eloglik_GP = function(x, s, Et, Et2) {
+  # Deal with infinite SEs:
+  idx = is.finite(s)
+  x = x[idx]
+  s = s[idx]
+  Et = Et[idx]
+  Et2 = Et2[idx]
+  return(-0.5 * sum(log(2*pi*s^2) + (1/s^2) * (Et2 - 2*x*Et + x^2)))
+}
+
+
+