Idea

The idea is to investigate method of moments estimation for very fast EB procedures like EBayesThresh.

Suppose \(x_j = b_j + N(0,s_j^2)\) and \(b_j \sim \pi_0 \delta_0 + (1-\pi_0) DExp(a)\) where \(Dexp\) means the double exponential (Laplace) with mean \(a\) (so rate parameter \(1/a\)).

Then the even moments of \(b_j\) are \(E(b_j^n) = (1-\pi_0) n! a^n\).

Also, 2nd and 4th moments of \(N(0,s_j^2)\) are \(s_j^2\) and \(3s_j^4\).

Thus, second and 4th moments of \(x_j\) are: \[E(x_j^2) = s_j^2 + 2(1-\pi_0)a^2\]

\[E(x_j^4) = 3s_j^4 + 24(1-\pi_0)a^4 + 6 s_j^2 [2(1-\pi_0)a^2] \].

Let \(m_2 := E(x_j^2 - s_j^2)\) and \(m_4:= E(x_j^4 - 3s_j^4)\). Then

\[E(m_2) = 2(1-\pi_0)a^2\]. \[E(m_4) = 24[(1-\pi_0)a^2] + 12s_j^2 (1-\pi_0)a^2\] \[= m_2 (12a^2 + 6s_j^2)\] So we can solve to give \[a^2 = (1/12) (m_4/m_2 - 6s_j^2)\]

OK, so I’ve been sloppy with subscripts on the \(m_4\) and \(m_2\)… need to sort that out.

We we will try it with constant \(s_j=1\).

Simulate some data

  eb_mm = function(x,s){
    m4 = mean(x^4-3*s^4)
    m2 = mean(x^2-s^2)
    a2 = ((m4/m2) - 6*mean(s^2))/12
    pi0 = 1- m2/(2*a2)
    if(a2<0){a2=0; pi0=1}
    if(pi0<0){pi0=0; a2 =m2/2}
    if(pi0>1){pi0=1; a2 = NA}
    
    return(list(pi0=pi0,a=sqrt(a2)))
  }
  set.seed(1)
  n=100000
  e = rnorm(n)
  b = rexp(n)
  x = b+e   
  eb_mm(x,1)

$pi0
[1] 0

$a
[1] 0.9969661

  EbayesThresh::wandafromx(x,1,FALSE)

$w
[1] 1

$a
[1] 1.003061

  b2 =b
  b2[1:(n/2)] = 0
  x2 = b2 + e
  eb_mm(x2,1)

$pi0
[1] 0.4807352

$a
[1] 0.988471

  EbayesThresh::wandafromx(x2,1,FALSE)

$w
[1] 0.5132097

$a
[1] 1.006022

Try a case that is “nearly null”

  b3 =b
  b3[1:(0.99*n)] = 0
  x3 = b3 + e
  eb_mm(x3,1)

$pi0
[1] 0.9822534

$a
[1] 0.8542569

  EbayesThresh::wandafromx(x3,1,FALSE)

$w
[1] 0.01192894

$a
[1] 1.062469

Thoughts

• hard part is when pi0 is very close to 1 but not 1. That might be worth thinking about. Maybe some upper quantiles would be better than 4th moment?

• the case where s is not constant - probably want to work with E(x/s) rather than E(x) ?

sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X El Capitan 10.11.6

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     

loaded via a namespace (and not attached):
 [1] MASS_7.3-47         backports_1.1.1     magrittr_1.5       
 [4] rprojroot_1.2       tools_3.3.2         htmltools_0.3.6    
 [7] yaml_2.1.14         Rcpp_0.12.13        stringi_1.1.5      
[10] rmarkdown_1.7       EbayesThresh_1.4-12 knitr_1.17         
[13] wavethresh_4.6.8    git2r_0.19.0        stringr_1.2.0      
[16] digest_0.6.12       evaluate_0.10.1