diff --git a/R/flash_data.R b/R/flash_data.R index 34f963e..e4f9340 100644 --- a/R/flash_data.R +++ b/R/flash_data.R @@ -9,7 +9,11 @@ #' ii) data$Y * data$missing is 0 if the original data were missing #' @return a flash data object #' @export -set_flash_data = function(Y, init=c("softimpute","mean")){ +set_flash_data = function(Y, init=c("softimpute","mean"), + noise_type = c("constant","known","Bayes_var","var_col","noisy","noisy_col"), + sigmae2_tilde = NA){ + noise_type = match.arg(noise_type, c("constant","known","Bayes_var","var_col","noisy","noisy_col")) + init = match.arg(init) data = list(Y=NULL, anyNA=FALSE, missing = FALSE) # initialize data @@ -25,6 +29,8 @@ set_flash_data = function(Y, init=c("softimpute","mean")){ Y[data$missing] = 0 } data$Y=Y + data$sigmae2_tilde = sigmae2_tilde + data$noise_type = noise_type return(data) } diff --git a/R/tau_est.R b/R/tau_est.R new file mode 100644 index 0000000..8f28efe --- /dev/null +++ b/R/tau_est.R @@ -0,0 +1,186 @@ + +#' title rescale the lambda_l and lambda_f for identifiablity +#' +#' description rescale the lambda_l and lambda_f for identifiablity +#' +#' @return a list of sig2_l and sig2_f for the rescaled variance of the kronecker product +#' @param sig2_l variance of the kronecker product +#' @param sig2_l variance of the kronecker product +#' @keywords internal +#' +rescale_sigmae2 = function(sig2_l,sig2_f){ + norm_l = sqrt(sum(sig2_l^2)) + norm_f = sqrt(sum(sig2_f^2)) + norm_total = norm_l * norm_f + sig2_l = sig2_l / norm_l + sig2_f = sig2_f / norm_f + sig2_l = sig2_l * sqrt(norm_total) + sig2_f = sig2_f * sqrt(norm_total) + return(list(sig2_l = sig2_l,sig2_f = sig2_f)) +} + +#' title initial value for Bayes variance structure estimation for kronecker productor +#' +#' description initial value for Bayes variance structure estimation for kronecker productor +#' +#' @return sig2_out estimated variance structure +#' @param sigmae2_v residual matrix +#' @keywords internal +#' +inital_Bayes_var = function(sigmae2_v){ + N = dim(sigmae2_v)[1] + P = dim(sigmae2_v)[2] + # estimate the initial value of lambda_l and lambda_f + sig2_l_pre = rowMeans(sigmae2_v, na.rm = TRUE) + sig2_f_pre = colMeans( sigmae2_v / matrix(rep(sig2_l_pre,P), ncol = P) , na.rm = TRUE) + sig2_pre_list = rescale_sigmae2(sig2_l_pre,sig2_f_pre) + sig2_l_pre = sig2_pre_list$sig2_l + sig2_f_pre = sig2_pre_list$sig2_f + # start the iteration + maxiter = 100 + inital_tol = 1e-3 + tau = 0 + epsilon = 1 + while(epsilon > inital_tol & tau <= maxiter){ + tau = tau + 1 + sig2_l = rowMeans( sigmae2_v / matrix(rep(sig2_f_pre,each = N),ncol = P) , na.rm = TRUE) + sig2_f = colMeans( sigmae2_v / matrix(rep(sig2_l,P), ncol = P) , na.rm = TRUE) + sig2_list = rescale_sigmae2(sig2_l,sig2_f) + sig2_l = sig2_list$sig2_l + sig2_f = sig2_list$sig2_f + epsilon = sqrt(mean((sig2_f - sig2_f_pre)^2 )) + sqrt(mean((sig2_l - sig2_l_pre)^2)) + sig2_l_pre = sig2_l + sig2_f_pre = sig2_f + } + return(list(sig2_l = sig2_l,sig2_f = sig2_f)) +} + +#' title Bayes variance structure estimation for kronecker productor +#' +#' description prior and posterior part in objective function +#' +#' @return sig2_out estimated variance structure +#' @param sigmae2_v residual matrix +#' @param sigmae2 variance structure +#' and it is a list of two vectors in this case. +#' @keywords internal +#' +Bayes_var = function(sigmae2_v,sigmae2 = NULL){ + N = dim(sigmae2_v)[1] + P = dim(sigmae2_v)[2] + if( is.null(sigmae2) || !is.list(sigmae2) ){ + # we don't know the truth this is in the first iteration + sigmae2 = inital_Bayes_var(sigmae2_v) + } + sig2_l_pre = sigmae2$sig2_l + sig2_f_pre = sigmae2$sig2_f + # this has already beed rescaled + # here we use alpha_l = alpha_f = beta_l = beta_f = 0 + sig2_l = rowMeans( sigmae2_v / matrix(rep(sig2_f_pre,each = N),ncol = P) , na.rm = TRUE) + sig2_f = colMeans( sigmae2_v / matrix(rep(sig2_l,P), ncol = P) , na.rm = TRUE) + #rescaled the variance + sig2_list = rescale_sigmae2(sig2_l,sig2_f) + sig2_l = sig2_list$sig2_l + sig2_f = sig2_list$sig2_f + # sig2_out = matrix(rep(sig2_l,P),ncol = P) * matrix(rep(sig2_f,each = N),ncol = P) + return(list(sig2_l = sig2_l,sig2_f = sig2_f)) +} + + +#' title noisy variance structure estimation +#' +#' description prior and posterior part in objective function +#' +#' @return sig2_out estimated variance structure +#' @param sigmae2_v residual matrix +#' @param sigmae2_true variance structure +#' and it is a list of two vectors in this case. +#' @keywords internal +#' + +noisy_var = function(sigmae2_v,sigmae2_true){ + # a reaonable upper bound which control optimal algorithm + upper_range = abs(mean(sigmae2_v - sigmae2_true , na.rm = TRUE)) + sd(sigmae2_v - sigmae2_true, na.rm = TRUE) + # value of likelihood + # I need the negative of the f_lik since optim find the minimum rather than the maximum + f_lik <- function(x,sigmae2_v_l = sigmae2_v,sigmae2_true_l = sigmae2_true){ + (1/2)*mean( log(2*pi*(x+sigmae2_true_l)) + (1/(x+sigmae2_true_l))*sigmae2_v_l , na.rm = TRUE) + } + AA <- optim(mean(sigmae2_v - sigmae2_true, na.rm = TRUE), f_lik,lower = 0, upper = upper_range, method = "Brent") + e_sig = AA$par + return(e_sig + sigmae2_true) +} + + +#' title noisy variance structure estimation with noisy variance on each column +#' +#' description prior and posterior part in objective function +#' +#' @return sig2_out estimated variance structure +#' @param sigmae2_v residual matrix +#' @param sigmae2_true variance structure +#' and it is a list of two vectors in this case. +#' @keywords internal +#' + +noisy_var_column = function(sigmae2_v,sigmae2_true){ + # in this case, we need to think about the for each column + # we can use noisy_var = function(sigmae2_v,sigmae2_true) with the input of each column + P = dim(sigmae2_true)[2] + sig2_out = sapply(seq(1,P),function(x){noisy_var(sigmae2_v[,x],sigmae2_true[,x])}) + return(sig2_out) +} + + +#' title module for estiamtion of the variance structure +#' +#' description estiamtion of the variance structure +#' +#' @return sigmae2 estimated variance structure +#' @param partype parameter type for the variance, +#' "constant" for constant variance, +#' "var_col" for nonconstant variance for column, +#' "known" for the kown variance, +#' "Bayes_var" for Bayes version of the nonconstant variance for row and column +#' "noisy" is noisy version s_ij + sigmae2 +#' "noisy_col" noisy version for column s_ij + sigmae2_j +#' @param sigmae2_v residual matrix +#' @param sigmae2_true true variance structure +#' @param sigmae2_pre the previouse sigmae2 which is only used in Bayes_var +#' @keywords internal +#' +# sigma estimation function +tau_est = function(f,data){ + # to get the residual matrix + sigmae2_v = get_R2(data,f) + sigmae2_v[data$missing] = NA + N = get_n(f) + P = get_p(f) + if(data$noise_type == "var_col"){ + # this case, we condier the sigma_j column-wise variance + sigmae2 = colMeans(sigmae2_v,na.rm = TRUE) + f$tau = 1/as.matrix(rep(1,N) %*% t(sigmae2)) + } else if(data$noise_type == "noisy"){ + # in this case we consider the sigmae2tilde + sigmae2 + sigmae2 = noisy_var(sigmae2_v,data$sigmae2_tilde) + f$tau = 1/sigmae2 + } else if(data$noise_type == "noisy_col"){ + # in this case we consider the sigmae2tilde + sigmae2_j + sigmae2 = noisy_var_column(sigmae2_v,data$sigmae2_tilde) + f$tau = 1/sigmae2 + }else if (data$noise_type == "Bayes_var"){ + # here sigmae2 is a list of two vectors and the variance is kronecker product of them + f$sigmae2 = Bayes_var(sigmae2_v,f$sigmae2) + sigmae2 = f$sigmae2$sig2_l %*% t(f$sigmae2$sig2_f) + f$tau = 1/sigmae2 + }else if (data$noise_type == "known"){ + sigmae2 = data$sigmae2_tilde + f$tau = 1/sigmae2 + } else { + # this is for constant case + sigmae2 = mean(sigmae2_v,na.rm = TRUE) + f$tau = 1/matrix(sigmae2,ncol = P, nrow = N) + } + return(f) +} + diff --git a/R/updates.R b/R/updates.R index 5106dd0..a0d7a82 100644 --- a/R/updates.R +++ b/R/updates.R @@ -67,7 +67,8 @@ flash_update_single_factor = function(data,f,k,ash_param=list()){ #' @title Update a single flash factor-loading combination (and precision) #' @inheritParams flash_update_single_loading flash_update_single_fl = function(data,f,k,ash_param=list()){ - f = flash_update_precision(data,f) + # f = flash_update_precision(data,f) + f = tau_est(f,data) f = flash_update_single_factor(data,f,k,ash_param) f = flash_update_single_loading(data,f,k,ash_param) return(f) @@ -132,5 +133,3 @@ flash_update_precision = function(data,f){ return(f) } - - diff --git a/man/Bayes_var.Rd b/man/Bayes_var.Rd new file mode 100644 index 0000000..6af0813 --- /dev/null +++ b/man/Bayes_var.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tau_est.R +\name{Bayes_var} +\alias{Bayes_var} +\title{title Bayes variance structure estimation for kronecker productor} +\usage{ +Bayes_var(sigmae2_v, sigmae2 = NULL) +} +\arguments{ +\item{sigmae2_v}{residual matrix} + +\item{sigmae2}{variance structure +and it is a list of two vectors in this case.} +} +\value{ +sig2_out estimated variance structure +} +\description{ +description prior and posterior part in objective function +} +\keyword{internal} diff --git a/man/inital_Bayes_var.Rd b/man/inital_Bayes_var.Rd new file mode 100644 index 0000000..50b456b --- /dev/null +++ b/man/inital_Bayes_var.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tau_est.R +\name{inital_Bayes_var} +\alias{inital_Bayes_var} +\title{title initial value for Bayes variance structure estimation for kronecker productor} +\usage{ +inital_Bayes_var(sigmae2_v) +} +\arguments{ +\item{sigmae2_v}{residual matrix} +} +\value{ +sig2_out estimated variance structure +} +\description{ +description initial value for Bayes variance structure estimation for kronecker productor +} +\keyword{internal} diff --git a/man/noisy_var.Rd b/man/noisy_var.Rd new file mode 100644 index 0000000..37e5076 --- /dev/null +++ b/man/noisy_var.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tau_est.R +\name{noisy_var} +\alias{noisy_var} +\title{title noisy variance structure estimation} +\usage{ +noisy_var(sigmae2_v, sigmae2_true) +} +\arguments{ +\item{sigmae2_v}{residual matrix} + +\item{sigmae2_true}{variance structure +and it is a list of two vectors in this case.} +} +\value{ +sig2_out estimated variance structure +} +\description{ +description prior and posterior part in objective function +} +\keyword{internal} diff --git a/man/noisy_var_column.Rd b/man/noisy_var_column.Rd new file mode 100644 index 0000000..3f73cff --- /dev/null +++ b/man/noisy_var_column.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tau_est.R +\name{noisy_var_column} +\alias{noisy_var_column} +\title{title noisy variance structure estimation with noisy variance on each column} +\usage{ +noisy_var_column(sigmae2_v, sigmae2_true) +} +\arguments{ +\item{sigmae2_v}{residual matrix} + +\item{sigmae2_true}{variance structure +and it is a list of two vectors in this case.} +} +\value{ +sig2_out estimated variance structure +} +\description{ +description prior and posterior part in objective function +} +\keyword{internal} diff --git a/man/rescale_sigmae2.Rd b/man/rescale_sigmae2.Rd new file mode 100644 index 0000000..c686f69 --- /dev/null +++ b/man/rescale_sigmae2.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tau_est.R +\name{rescale_sigmae2} +\alias{rescale_sigmae2} +\title{title rescale the lambda_l and lambda_f for identifiablity} +\usage{ +rescale_sigmae2(sig2_l, sig2_f) +} +\arguments{ +\item{sig2_l}{variance of the kronecker product} + +\item{sig2_l}{variance of the kronecker product} +} +\value{ +a list of sig2_l and sig2_f for the rescaled variance of the kronecker product +} +\description{ +description rescale the lambda_l and lambda_f for identifiablity +} +\keyword{internal} diff --git a/man/set_flash_data.Rd b/man/set_flash_data.Rd index 96bc46a..7a8e80d 100644 --- a/man/set_flash_data.Rd +++ b/man/set_flash_data.Rd @@ -4,7 +4,8 @@ \alias{set_flash_data} \title{set up data for reading into flash} \usage{ -set_flash_data(Y, init = c("softimpute", "mean")) +set_flash_data(Y, init = c("softimpute", "mean"), noise_type = c("constant", + "known", "Bayes_var", "var_col", "noisy", "noisy_col"), sigmae2_tilde = NA) } \arguments{ \item{Y}{an n by p data matrix} diff --git a/man/tau_est.Rd b/man/tau_est.Rd new file mode 100644 index 0000000..f0e4a9b --- /dev/null +++ b/man/tau_est.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tau_est.R +\name{tau_est} +\alias{tau_est} +\title{title module for estiamtion of the variance structure} +\usage{ +tau_est(f, data) +} +\arguments{ +\item{partype}{parameter type for the variance, +"constant" for constant variance, +"var_col" for nonconstant variance for column, +"known" for the kown variance, +"Bayes_var" for Bayes version of the nonconstant variance for row and column +"noisy" is noisy version s_ij + sigmae2 +"noisy_col" noisy version for column s_ij + sigmae2_j} + +\item{sigmae2_v}{residual matrix} + +\item{sigmae2_true}{true variance structure} + +\item{sigmae2_pre}{the previouse sigmae2 which is only used in Bayes_var} +} +\value{ +sigmae2 estimated variance structure +} +\description{ +description estiamtion of the variance structure +} +\keyword{internal} diff --git a/vignettes/flash_r1_tau.Rmd b/vignettes/flash_r1_tau.Rmd new file mode 100644 index 0000000..4d37c69 --- /dev/null +++ b/vignettes/flash_r1_tau.Rmd @@ -0,0 +1,168 @@ +--- +title: "test the variance estimation" +author: "Wei Wang" +date: "9/19/2017" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +```{r} +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*(k/K)) + fstart = rnorm(P, 0, signal*(k/K)) + + 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)) +} + +sim_hd = function(N, P, SF, SL, signal, a = rchisq(N,3),b = rchisq(P,1),mu = 0){ + + 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(mu + a[i] * b[j])) + } + } + 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 + + Y = lstart %*% t(fstart) + E + + return(list(Y = Y, L_true = lstart, F_true = fstart, Error = E,sig2_true = sig2_true)) +} + +``` + + +```{r} +library(flashr2) +``` + +### constant variance + +```{r} +set.seed(99) +N = 100 +P = 200 +dat = sim_K(K=1,N, P, SF = 0.8, SL = 0.5, signal = 1,noise = 0.5) +Y = dat$Y +E = dat$Error + +data = set_flash_data(Y) +f= flash_r1(data,verbose=F) + +``` + +this version is already tested + +### non-constant variance + +In the following situations, we simulate the data with rank one structure with sparsity on both factor and loading. +The variance of error term is +$$\sigma_{e_{ij}}^2 = \mu + a_ib_j$$ + +#### columnwise variance + +we set $\mu = 0$ and $a = \bf{1}$, so the $\sigma_{e_{ij}}^2 = b_j$ + +```{r} +set.seed(99) +N = 100 +P = 200 +dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rep(1,N), b = rchisq(P,1), mu = 0) +Y = dat$Y +E = dat$Error + +data = set_flash_data(Y,noise_type = "var_col") +f= flash_r1(data,verbose=F) + +plot(dat$sig2_true[1,],as.vector(1/f$tau[1,])) +abline(0,1) +``` + +#### Kronecker product variance + +we set $\mu = 0$, so the $\sigma_{e_{ij}}^2 = a_ib_j$ + +```{r} +set.seed(99) +N = 100 +P = 200 +dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rchisq(N,3),b = rchisq(P,1),mu = 0) +Y = dat$Y +E = dat$Error + +data = set_flash_data(Y,noise_type = "B") +f= flash_r1(data,verbose=F) + + +plot(as.vector(dat$sig2_true),as.vector(1/f$tau)) +abline(0,1) +``` + +#### Known variance matrix + +```{r} +set.seed(99) +N = 100 +P = 200 +dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rchisq(N,3),b = rchisq(P,1),mu = 0) +Y = dat$Y +E = dat$Error + +data = set_flash_data(Y,noise_type = "k",sigmae2_tilde = dat$sig2_true) +f= flash_r1(data,verbose=F) + + +plot(as.vector(dat$sig2_true),as.vector(1/f$tau)) +abline(0,1) +``` + + +#### Noisy variance + +we set $\mu = 0.9$, so the $\sigma_{e_{ij}}^2 = a_ib_j + 0.9$. We assume that $a_ib_j$ is known, so we set `sigmae2_tilde = (dat$sig2_true - 0.9)` + +To test the estimation, we use `1/f$tau - sigmae2_tilde` where `sigmae2_tilde = (dat$sig2_true - 0.9) = (a_ib_j)`. We expect `1/f$tau - sigmae2_tilde` to be a matrix of a constant around 0.9. + + +```{r} +set.seed(99) +N = 100 +P = 200 +dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rchisq(N,3),b = rchisq(P,1),mu = 0.9) +Y = dat$Y +E = dat$Error + +data = set_flash_data(Y,noise_type = "noisy",sigmae2_tilde = (dat$sig2_true - 0.9)) +f= flash_r1(data,verbose=F) + +(1/f$tau - (dat$sig2_true - 0.9))[1:10,1:10] + +``` + diff --git a/vignettes/flash_r1_tau.html b/vignettes/flash_r1_tau.html new file mode 100644 index 0000000..fd637b4 --- /dev/null +++ b/vignettes/flash_r1_tau.html @@ -0,0 +1,309 @@ + + + + + + + + + + + + + + +test the variance estimation + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + +
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*(k/K))
+    fstart = rnorm(P, 0, signal*(k/K))
+    
+    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))
+}
+
+sim_hd = function(N, P, SF, SL, signal, a = rchisq(N,3),b = rchisq(P,1),mu = 0){
+  
+  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(mu + a[i] * b[j]))
+    }
+  }
+  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
+  
+  Y = lstart %*% t(fstart) + E
+  
+  return(list(Y = Y, L_true = lstart, F_true = fstart, Error = E,sig2_true = sig2_true))
+}
+
library(flashr2)
+
+

constant variance

+
set.seed(99)
+N = 100
+P = 200
+dat = sim_K(K=1,N, P, SF = 0.8, SL = 0.5, signal = 1,noise = 0.5)
+Y = dat$Y
+E = dat$Error
+
+data = set_flash_data(Y)
+f= flash_r1(data,verbose=F)
+

this version is already tested

+
+
+

non-constant variance

+

In the following situations, we simulate the data with rank one structure with sparsity on both factor and loading. The variance of error term is \[\sigma_{e_{ij}}^2 = \mu + a_ib_j\]

+
+

columnwise variance

+

we set \(\mu = 0\) and \(a = \bf{1}\), so the \(\sigma_{e_{ij}}^2 = b_j\)

+
set.seed(99)
+N = 100
+P = 200
+dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rep(1,N), b = rchisq(P,1), mu = 0)
+Y = dat$Y
+E = dat$Error
+
+data = set_flash_data(Y,noise_type = "var_col")
+f= flash_r1(data,verbose=F)
+
+plot(dat$sig2_true[1,],as.vector(1/f$tau[1,]))
+abline(0,1)
+

+
+
+

Kronecker product variance

+

we set \(\mu = 0\), so the \(\sigma_{e_{ij}}^2 = a_ib_j\)

+
set.seed(99)
+N = 100
+P = 200
+dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rchisq(N,3),b = rchisq(P,1),mu = 0)
+Y = dat$Y
+E = dat$Error
+
+data = set_flash_data(Y,noise_type = "B")
+f= flash_r1(data,verbose=F)
+
+
+plot(as.vector(dat$sig2_true),as.vector(1/f$tau))
+abline(0,1)
+

+
+
+

Known variance matrix

+
set.seed(99)
+N = 100
+P = 200
+dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rchisq(N,3),b = rchisq(P,1),mu = 0)
+Y = dat$Y
+E = dat$Error
+
+data = set_flash_data(Y,noise_type = "k",sigmae2_tilde =  dat$sig2_true)
+f= flash_r1(data,verbose=F)
+
+
+plot(as.vector(dat$sig2_true),as.vector(1/f$tau))
+abline(0,1)
+

+
+
+

Noisy variance

+

we set \(\mu = 0.9\), so the \(\sigma_{e_{ij}}^2 = a_ib_j + 0.9\). We assume that \(a_ib_j\) is known, so we set sigmae2_tilde = (dat$sig2_true - 0.9)

+

To test the estimation, we use 1/f$tau - sigmae2_tilde where sigmae2_tilde = (dat$sig2_true - 0.9) = (a_ib_j). We expect 1/f$tau - sigmae2_tilde to be a matrix of a constant around 0.9.

+
set.seed(99)
+N = 100
+P = 200
+dat = sim_hd(N, P, SF = 0.5, SL=0.5, signal =1, a = rchisq(N,3),b = rchisq(P,1),mu = 0.9)
+Y = dat$Y
+E = dat$Error
+
+data = set_flash_data(Y,noise_type = "noisy",sigmae2_tilde =  (dat$sig2_true - 0.9))
+f= flash_r1(data,verbose=F)
+
+(1/f$tau - (dat$sig2_true - 0.9))[1:10,1:10]
+
##            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
+##  [1,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [2,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [3,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [4,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [5,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [6,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [7,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [8,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##  [9,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+## [10,] 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546 0.8971546
+##            [,7]      [,8]      [,9]     [,10]
+##  [1,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [2,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [3,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [4,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [5,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [6,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [7,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [8,] 0.8971546 0.8971546 0.8971546 0.8971546
+##  [9,] 0.8971546 0.8971546 0.8971546 0.8971546
+## [10,] 0.8971546 0.8971546 0.8971546 0.8971546
+
+
+ + + + +
+ + + + + + + +