Last updated: 2017-06-20

Code version: 177286a

In this case, we would like to try the objective function with penalty term.

I would like to modify the flash code in this document.

sim_K = function(K, N, P, SF, SL, signal,noise){
  E = matrix(rnorm(N*P,0,noise),nrow=N)
  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
    
    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))
}

check the new obj function

we check that the new objective function with penalty term does monotonically increase.

set.seed(99)
objcheck_1 = function(K=1,N = 20, P=30, SF = 0.5, SL = 0.6, signal = 1,noise = 1){
  par(mfrow = c(3, 2))
  par(cex = 0.6)
  par(mar = c(3, 3, 0.8, 0.8), oma = c(1, 1, 1, 1))
  data = sim_K(K,N, P, SF , SL , signal ,noise )
  Y = data$Y
  E = data$Error
  g1 = flash(Y,objtype = "l")
  plot(g1$obj_val_track,main = "fdr method with penalized obj",type = "l")
  g2 =  flash(Y,objtype = "l",ash_para = list(method = "shrink"))
  plot(g2$obj_val_track,main = "shrink method without penalized obj",type = "l")
  g3 =  flashr::flash(Y,objtype = "l",ash_para = list(method = "shrink"))
  plot(g3$obj_val_track,main = "shrink method without penalized obj(old)",type = "l")
  g4 =  flashr::flash(Y,objtype = "l")
  plot(g4$obj_val_track,main = "fdr method without penalized obj",type = "l")
  plot(c(g1$obj_val,g2$obj_val,g3$obj_val,g4$obj_val),main = "obj value")
  RMSE = c(sqrt(mean((Y - g1$l %*% t(g1$f) -E)^2)) ,
           sqrt(mean((Y - g2$l %*% t(g2$f) -E)^2)) ,
           sqrt(mean((Y - g3$l %*% t(g3$f) -E)^2)) ,
           sqrt(mean((Y - g4$l %*% t(g4$f) -E)^2)) )
  plot(RMSE,main = "RMSE")
}
objcheck_1(K=1,N = 20, P=30, SF = 0.5, SL = 0.6, signal = 1,noise = 1)

objcheck_1(K=1,N = 200, P=300, SF = 0.8, SL = 0.6, signal = 1,noise = 1)

objcheck_1(K=1,N = 200, P=300, SF = 0.9, SL = 0.8, signal = 1,noise = 1)

objcheck_1(K=1,N = 200, P=300, SF = 0.5, SL = 0.9, signal = 1,noise = 3)

objcheck_1(K=1,N = 100, P=200, SF = 0.8, SL = 0.9, signal = 1,noise = 3)

we can see that the old version of code without penalty doesn’t have monotonically increasing objective function, but the new version code with penalty term does.

now we have two methods (fdr, shrink), and both of them have monotonically increasing corresponding objective function in each case.

try the same problem: rank 0 vs rank 1 with fdr method.

case 1: rank 0

we compare the old code with our new code, we can see that most of the cases a

shrink method without penalty

set.seed(99)
sim_rank0(K=1,N = 100, P=200, SF = 0.5, SL = 0.5, signal = 0,noise = 1,mtype = "shrink")

[1] 0.400 0.005 0.595

fdr method with penalty

set.seed(99)
sim_rank0(K=1,N = 100, P=200, SF = 0.5, SL = 0.5, signal = 0,noise = 1)

[1] 0.875 0.000 0.125

we can see the improvement that the proportion of P(estimated structure better than rank 0 structure) is much larger.

case 2: rank 1, very sparse with large noise

shrink method without penalty

set.seed(99)
sim_rank0(K=1,N = 100, P=200, SF = 0.8, SL = 0.8, signal = 1,noise = 2,mtype = "shrink")

[1] 0.100 0.455 0.445

fdr method with penalty

set.seed(99)
sim_rank0(K=1,N = 100, P=200, SF = 0.8, SL = 0.8, signal = 1,noise = 2)

[1] 0.280 0.415 0.305

we can see that the results are improved since the proportion of P(estimated structure better than rank 0 structure) is much larger.

In this case we can see that the RMSE is negative correlated with the objective function value. So it is make sense that we return \(g_l = 0, g_f = 0\) if the \(LB(estimated) < LB(g_l=0,g_f = 0)\).

Session information

sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.12.4 (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] workflowr_0.4.0 rmarkdown_1.3  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.11      knitr_1.15.1      magrittr_1.5     
 [4] REBayes_0.73      MASS_7.3-45       munsell_0.4.3    
 [7] doParallel_1.0.10 pscl_1.4.9        colorspace_1.3-2 
[10] SQUAREM_2016.8-2  lattice_0.20-34   foreach_1.4.3    
[13] plyr_1.8.4        flashr_0.1.1      ashr_2.1-21      
[16] stringr_1.2.0     tools_3.3.0       parallel_3.3.0   
[19] grid_3.3.0        gtable_0.2.0      irlba_2.1.2      
[22] git2r_0.18.0      htmltools_0.3.5   iterators_1.0.8  
[25] lazyeval_0.2.0    assertthat_0.2.0  yaml_2.1.14      
[28] rprojroot_1.2     digest_0.6.12     tibble_1.2       
[31] Matrix_1.2-8      ggplot2_2.2.1     codetools_0.2-15 
[34] evaluate_0.10     stringi_1.1.2     scales_0.4.1     
[37] Rmosek_7.1.2      backports_1.0.5   truncnorm_1.0-7

This R Markdown site was created with workflowr