Processing math: 100%
  • Idea
  • Simulate some data
  • Thoughts
    • Session information

Last updated: 2017-12-12

Code version: ce26fc9

Idea

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

Suppose xj=bj+N(0,s2j) and bjπ0δ0+(1π0)DExp(a) where Dexp means the double exponential (Laplace) with rate a.

Then the even moments of bj are E(bnj)=(1π0)n!(1/a)n.

Also, 2nd and 4th moments of N(0,s2j) are s2j and 3s4j.

Thus, second and 4th moments of xj are: E(x2j)=s2j+2(1π0)/a2

E(x4j)=3s4j+24(1π0)a4+6s2j[2(1π0)a2].

Let m2:=E(x2js2j) and m4:=E(x4j3s4j). Then

E(m2)=2(1π0)/a2. E(m4)=24[(1π0)/a2]+12s2j(1π0)/a2 =m2(12/a2+6s2j) So we can solve to give a2=12/(m4/m26s2j)

OK, so I’ve been sloppy with subscripts on the m4 and m2… need to sort that out.

We we will try it with constant sj=1.

Simulate some data

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

$a
[1] 2.003444
  EbayesThresh::wandafromx(x,1,FALSE)
$w
[1] 1

$a
[1] 2.00342
  b2 =b
  b2[1:(n/2)] = 0
  x2 = b2 + e
  eb_mm(x2,1)
$pi0
[1] 0.4384081

$a
[1] 2.086427
  EbayesThresh::wandafromx(x2,1,FALSE)
$w
[1] 0.5581188

$a
[1] 2.079653
  EbayesThresh:::wandafromx.mle(x2,1)
$w
[1] 0.558119

$a
[1] 2.079654

Try a case that is “nearly null”. Note that here the original Ebayesthresh approach based on the beta function is less accurate, presumably due to numeric issues.

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

$a
[1] 3.961431
  EbayesThresh::wandafromx(x3,1,FALSE)
$w
[1] 1

$a
[1] 13.08551
  EbayesThresh:::wandafromx.mle(x3,1)
$w
[1] 0.09129917

$a
[1] 3.961397

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) ?

Session information

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.14       stringi_1.1.5     
[10] rmarkdown_1.7      EbayesThresh_1.5-0 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   

This R Markdown site was created with workflowr