Last updated: 2016-10-25

Code version: 3d29620508ebbe9aa0d06740c6588ac415e9ca1d

generate data

sim_data = function(N, P, SF, SL, signal,
                    a = rchisq(N,3), b = rchisq(P,1),
                    mu = 1, K = 6, positive = FALSE){
  
  E = matrix(rep(0,N*P),nrow=N)
  sig2_true = matrix(rep(0,N*P),nrow=N)
  for(i in 1:N){
    for(j in 1:P){
      sig2_true[i,j] = mu + a[i] * b[j]
      E[i,j] = rnorm(1,0,sqrt(sig2_true[i,j]))
    }
  }
  
  Y = E
  L_true = array(0, dim = c(N,K))
  F_true = array(0, dim = c(P,K))
  
  for(k in 1:K){
    lstart = rnorm(N, 0, signal)
    fstart = rnorm(P, 0, signal)
    
    index = sample(seq(1:N),(N*SL))
    lstart[index] = 0
    index = sample(seq(1:P),(P*SF))
    fstart[index] = 0
    
    if(positive == TRUE){
      lstart = abs(lstart)
      fstart = abs(fstart)
    }
    
    L_true[,k] = lstart
    F_true[,k] = fstart
    
    Y = Y + lstart %*% t(fstart)
  }
  
  return(list(Y = Y, L_true = L_true, F_true = F_true, Error = E,sig2_true = sig2_true))
}

try the rank one case

we set the variance is a constant plus a columnwise structure variance.

set.seed(100)
N = 60
P = 100
SF = 0.5
SL = 0.8
data = sim_data(N, P, SF, SL, a = rep(1,N), b = rchisq(P,1), mu = 1.5, signal = 1,K = 1)

plot(svd(data$Y)$d)
# data$sig2_true

# First test the flash rank one case.
g_flash = flashr::flash(data$Y,partype = "noisy_col",sigmae2_true = matrix(1.5,ncol = P, nrow = N))

plot((data$sig2_true[1,]), g_flash$sigmae2[1,])
abline(0,1)

set.seed(100)
N = 60
P = 100
SF = 0.5
SL = 0.8
data = sim_data(N, P, SF, SL,a = rep(1,N), b = rchisq(P,1), mu = 0.2, signal = 2,K = 1)

plot(svd(data$Y)$d)

# data$sig2_true

# First test the flash rank one case.
g_flash = flashr::flash(data$Y,partype = "noisy_col",sigmae2_true = matrix(0.2,ncol = P, nrow = N))
plot((data$sig2_true[1,]), g_flash$sigmae2[1,])
abline(0,1)

known constant plus a unknown constant

set.seed(100)
N = 60
P = 100
SF = 0.5
SL = 0.8
data = sim_data(N, P, SF, SL,a = rep(1,N), b = rep(2,P), mu = 0.2, signal = 2,K = 1)

plot(svd(data$Y)$d)

g_flash = flashr::flash(data$Y,partype = "noisy",sigmae2_true = matrix(0.2,ncol = P, nrow = N))
c((data$sig2_true[1,1]), g_flash$sigmae2[1,1])
[1] 2.200000 2.178848

test for the rank K, for the backfitting algorithm we only use 10 iterations

known constant plus unknown constant

set.seed(100)
N = 60
P = 100
SF = 0.5
SL = 0.8
data = sim_data(N, P, SF, SL,a = rep(1,N), b = rep(2,P), mu = 1, signal = 2,K = 6)
plot(svd(data$Y)$d)

g_flash = flashr::greedy(data$Y,K = 9,
                         flash_para = list(partype = "noisy",sigmae2_true = matrix(1,ncol = P, nrow = N)))
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
b_flash = flashr::backfitting(data$Y,initial_list = g_flash, maxiter_bf = 10, 
                              flash_para = list(partype = "noisy",sigmae2_true = matrix(1,ncol = P, nrow = N)))

c(data$sig2_true[1,1],(b_flash$sigmae2[[6]])[1,1])
[1] 3.000000 2.981137

known constant plus unknows columnwise structure

set.seed(100)
N = 60
P = 100
SF = 0.5
SL = 0.8
data = sim_data(N, P, SF, SL,a = rep(1,N), b = rchisq(P,1), mu = 1, signal = 2,K = 6)
plot(svd(data$Y)$d)

g_flash = flashr::greedy(data$Y,K = 9,
                         flash_para = list(partype = "noisy_col",sigmae2_true = matrix(1,ncol = P, nrow = N)))
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
b_flash = flashr::backfitting(data$Y,initial_list = g_flash, maxiter_bf = 10, 
                              flash_para = list(partype = "noisy_col",sigmae2_true = matrix(1,ncol = P, nrow = N)))
plot((data$sig2_true[1,]), (b_flash$sigmae2[[6]])[1,])
abline(0,1)

Session information

R version 3.2.4 (2016-03-10)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.12 (unknown)

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] knitr_1.14

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.7       magrittr_1.5      REBayes_0.63     
 [4] MASS_7.3-45       doParallel_1.0.10 pscl_1.4.9       
 [7] SQUAREM_2016.8-2  lattice_0.20-34   foreach_1.4.3    
[10] ashr_1.2.5        stringr_1.1.0     flashr_0.1.1     
[13] tools_3.2.4       parallel_3.2.4    grid_3.2.4       
[16] irlba_2.1.2       htmltools_0.3.5   iterators_1.0.8  
[19] assertthat_0.1    yaml_2.1.13       digest_0.6.10    
[22] Matrix_1.2-7.1    formatR_1.4       codetools_0.2-15 
[25] evaluate_0.10     rmarkdown_1.0     stringi_1.1.2    
[28] Rmosek_7.1.2      truncnorm_1.0-7