Result15_Estimation of Sigma

Introduction

This documentation is to analyze the effect of parametrization. Our convention is to write $y\sim g$ if $g$ is a probability density (with respect to a fixed measure). Let

$$g(\beta) = \sum_{k=1}^K \pi_k g_k(\beta)$$ is fixed and independent of $\sigma^2$. Let’s simply consider the normal case when $g_k = \mathcal{N}(0,\sigma_k^2)$.

1. Parametrization 1: This is what VARBVS uses. $$y|X,\beta,\sigma^2 \sim \mathcal{N}(X\beta,\sigma^2 I_n),\quad \beta_j | g, \sigma^2 = \sum_{k=1}^K \pi_k \mathcal{N}(0,\sigma^2\sigma_k^2), \quad q_j(\beta_j) = \sum_{k=1}^K \phi_{jk} \mathcal{N}(m_{jk},s_{jk}^2)$$ The update of $\beta_j$ is based on $$q_j^* = \arg\min_{q_j} \frac{(X^T X)_{jj}}{2\sigma^2}\mathbb{E}_{q_j} (\beta_j - b_j)^2 + D_{KL}(q_j\|g)$$ which is a posterior of $$b_j|\beta_j,\sigma^2 \sim \mathcal{N}(\beta_j,\sigma^2 / (X^T X)_{jj}),\quad \beta_j|g \sim \sum_k \pi_k \mathcal{N}(0,\sigma_k^2)$$ The update of $\sigma^2$ in this case is $$\sigma^2 = \frac{\|y - {\bf X}\mu\|^2 + \sum_j ({\bf X}^T{\bf X})_{jj} \mathrm{Var}[\beta_j] + \sum_j \sum_{k=2}^K \phi_{jk} (m_{jk}^2 + s_{jk}^2)/\sigma_k^2} {n + \sum_j \sum_{k=2}^K \phi_{jk}} = \frac{\|y - {\bf X}\mu\|^2 + \sum_j ({\bf X}^T{\bf X})_{jj} (\sum_{k} \phi_{jk} (m_{jk}^2 + s_{jk}^2) - \mu_j^2) + \sum_j \sum_{k=2}^K \phi_{jk} (m_{jk}^2 + s_{jk}^2)/\sigma_k^2} {n + \sum_j \sum_{k=2}^K \phi_{jk}}.$$ This is the varbvsmix’s update of $\sigma^2$. Compared to the update of Parametrization 2 below, it has additional $s_{jk}^2$ terms on the numerator and the $\sum_j \sum_{k=2}^K \phi_{jk}$ term on the denominator.

2. Parametrization 2: $$y|X,\beta,\sigma^2 \sim \mathcal{N}(X\beta,\sigma^2 I_n),\quad \beta_j | g, \sigma^2 = \sum_{k=1}^K \pi_k \mathcal{N}(0,\sigma^2\sigma_k^2), \quad q_j(\beta_j) = \sum_{k=1}^K \phi_{jk} \mathcal{N}(\mu_{jk},\sigma^2s_{jk}^2)$$ The update of $\beta_j$ is based on $$q_j^* = \arg\min_{q_j} \frac{(X^T X)_{jj}}{2\sigma^2}\mathbb{E}_{q_j} (\beta_j - b_j)^2 + D_{KL}(q_j\|g)$$ which is a posterior of $$b_j|\beta_j,\sigma^2 \sim \mathcal{N}(\beta_j,\sigma^2 / (X^T X)_{jj}),\quad \beta_j|g \sim \sum_k \pi_k \mathcal{N}(0,\sigma^2\sigma_k^2)$$ where $b_j = \mathbb{E}_q \beta_j + (X^T X)_{jj}^{-1}(y - X\beta)$. The update of $\sigma^2$ in this case is $$\sigma^2 = \frac{\|y - {\bf X}\mu\|^2 + \sum_j ({\bf X}^T{\bf X})_{jj} (\sum_{k} \phi_{jk} m_{jk}^2 - \mu_j^2) + \sum_j \sum_{k=2}^K \phi_{jk} m_{jk}^2/\sigma_k^2} {n}.$$ This can be reduced to $$\sigma^2 = \frac{\|y - {\bf X}\mu\|^2 + \sum_j ({\bf X}^T{\bf X})_{jj} (b_j\mu_j - \mu_j^2)} {n}.$$

3. Parametrization 3: $$y|X,\beta,\sigma^2 \sim \mathcal{N}(\sigma X\beta,\sigma^2 I_n),\quad \beta_j | g = \sum_{k=1}^K \pi_k \mathcal{N}(0,\sigma_k^2), \quad q_j(\beta_j) = \sum_{k=1}^K \phi_{jk} \mathcal{N}(\mu_{jk},s_{jk}^2)$$ The update of $\beta_j$ is based on $$q_j^* = \arg\min_{q_j} \frac{(X^T X)_{jj}}{2\sigma^2}\mathbb{E}_{q_j} (\beta_j - b_j)^2 + D_{KL}(q_j\|g)$$ which is a posterior of $$b_j|\beta_j,\sigma^2 \sim \mathcal{N}(\beta_j,1 / (X^T X)_{jj}),\quad \beta_j|g \sim \sum_k \pi_k \mathcal{N}(0,\sigma_k^2)$$ where $b_j = \mathbb{E}_q \beta_j + (X^T X)_{jj}^{-1}(y / \sigma - X\beta)$. In this case, the update of $\sigma^2$ is $$\sigma^2 = \arg\min_{\sigma} \frac{n}{2} \log(\sigma^2) + \frac{1}{2\sigma^2}\mathbb{E}_q \|y - \sigma X\beta \|^2 = \arg \min_\sigma \frac{n}{2} \log(\sigma^2) + \frac{1}{2\sigma^2} (y^T y - 2\sigma y^T X \mu) = \arg \min_\sigma A\log(\sigma^2) - \frac{2B}{\sigma} + \frac{C}{\sigma^2}$$ where $A = n$, $B = y^T X\mu$ and $C = y^T y$. Thus $\sigma$ is a solution to $$2A \sigma^2 + 2B \sigma - 2C = 0 \iff \sigma = \frac{-B + \sqrt{B^2 + 4AC}}{2A} = \frac{-(y^T X\mu) + \sqrt{(y^T X\mu)^2 + 4n(y^T y)}}{2n}$$

Results

tdat = readRDS("results/estsigma.RDS")
res = matrix(0,4,3)
for (i in 1:4) {
  res[i,] = colMeans(matrix(tdat[[i]]$sigma, 20, 3))
}
colnames(res) = c("original","Matthew","varbvsmix")
rownames(res) = c("pve = 0.1","pve = 0.5","pve = 0.9","pve = 0.9999")
print("average value of (sigma_est / sigma_true)")

[1] "average value of (sigma_est / sigma_true)"

res

             original   Matthew varbvsmix
pve = 0.1    1.061100 1.0610999  1.061073
pve = 0.5    1.039929 1.0399285  1.039890
pve = 0.9    1.003104 1.0030930  1.003073
pve = 0.9999 1.081021 0.7813146  3.597446

res = matrix(0,4,3)
for (i in 1:4) {
  res[i,] = colMeans(matrix(tdat[[i]]$rate, 20, 3))
}
colnames(res) = c("original","Matthew","varbvsmix")
rownames(res) = c("pve = 0.1","pve = 0.5","pve = 0.9","pve = 0.9999")
print("average value of min(sigma_est / sigma_true, sigma_true / sigma_est)")

[1] "average value of min(sigma_est / sigma_true, sigma_true / sigma_est)"

res

              original   Matthew varbvsmix
pve = 0.1    0.9269040 0.9269038 0.9269171
pve = 0.5    0.9318945 0.9318945 0.9318989
pve = 0.9    0.9291743 0.9291845 0.9292332
pve = 0.9999 0.9247187 0.7813146 0.3561272

res = matrix(0,4,3)
for (i in 1:4) {
  res[i,] = colMeans(matrix(tdat[[i]]$pred, 20, 3))
}
colnames(res) = c("original","Matthew","varbvsmix")
rownames(res) = c("pve = 0.1","pve = 0.5","pve = 0.9","pve = 0.9999")
print("prediction error (RMSE_test / sigma_true)")

[1] "prediction error (RMSE_test / sigma_true)"

res

             original    Matthew varbvsmix
pve = 0.1    1.040558   1.040652  1.040556
pve = 0.5    1.059430   1.060562  1.059428
pve = 0.9    1.035987   1.037942  1.036011
pve = 0.9999 1.012091 112.557985  1.014590

library(Matrix); library(ggplot2); library(cowplot); library(mr.ash); library(varbvs)
standardize = FALSE
source('code/method_wrapper.R')
source('code/sim_wrapper.R')
n   = 500
p   = 2000
s   = 20
sa2 = (2^((0:19)/5) - 1)^2
K   = length(sa2)
time = matrix(0,20,3)
sigma = matrix(0,20,3)
rate = matrix(0,20,3)
pred = matrix(0,20,3)
pve_list = c(0.1,0.5,0.9,0.9999)
tdat = list()

for (iter in 1:4) {
for (i in 1:20) {
  data          = simulate_data(n, p, s = s, seed = i, signal = "normal",
                                design = "indepgauss", pve = pve_list[iter])
  
  #sigma.init = var(data$y)
  sigma.init = data$sigma^2
  
  t.mr.ash1           = system.time(
    fit.mr.ash1        <- mr.ash(X = data$X, y = data$y, sa2 = sa2, sigma2 = sigma.init,
                                max.iter = 2000, min.iter = 1,
                                standardize = standardize,
                                tol = list(epstol = 1e-12, convtol = 1e-8)))
  
  t.mr.ash2           = system.time(
    fit.mr.ash2        <- mr.ash(X = data$X, y = data$y, sa2 = sa2, method = "em2",
                                max.iter = 2000, min.iter = 1, sigma2 = sigma.init,
                                standardize = standardize,
                                tol = list(epstol = 1e-12, convtol = 1e-8)))
  
  t.varbvsmix         = system.time(
    fit.varbvsmix      <- varbvsmix(X = data$X, Z = NULL, y = data$y, sa = sa2,
                                    mu = matrix(0, p, K), sigma = sigma.init,
                                    alpha = matrix(1, p, K) / K, update.sa = FALSE,
                                    update.sigma = TRUE, verbose = FALSE))
  
  pred1 = data$y.test - predict(fit.mr.ash1, data$X.test)
  pred2 = data$y.test - data$X.test %*% fit.mr.ash2$beta * sqrt(fit.mr.ash2$sigma) - fit.mr.ash2$intercept
  pred3 = data$y.test - data$X.test %*% c(rowSums(fit.varbvsmix$alpha * fit.varbvsmix$mu)) - fit.varbvsmix$mu.cov
  
  pred[i,] <- c(norm(pred1, '2'), norm(pred2, '2'), norm(pred3, '2')) / sqrt(n) / data$sigma
  a = fit.mr.ash1$sigma / data$sigma^2
  b = fit.mr.ash2$sigma / data$sigma^2
  c = fit.varbvsmix$sigma / data$sigma^2
  time[i,] = c(t.mr.ash1[3], t.mr.ash2[3], t.varbvsmix[3])
  rate[i,] = c("caisa" = min(a,1/a), "new" = min(b, 1/b), "varbvsmix" = min(c, 1/c))
  sigma[i,] = c(a,b,c)
}
tdat[[iter]] = data.frame(pred = c(pred), time = c(time), rate = c(rate), sigma = c(sigma),
                          fit = c("caisa","new","varbvsmix"))
print(colMeans(pred))
print(colMeans(rate))
}

Session information

sessionInfo()

R version 3.5.3 (2019-03-11)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

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] workflowr_1.4.0 Rcpp_1.0.2      digest_0.6.21   rprojroot_1.3-2
 [5] backports_1.1.4 git2r_0.26.1    magrittr_1.5    evaluate_0.14  
 [9] stringi_1.4.3   fs_1.3.1        rmarkdown_1.15  tools_3.5.3    
[13] stringr_1.4.0   glue_1.3.1      xfun_0.9        yaml_2.2.0     
[17] compiler_3.5.3  htmltools_0.3.6 knitr_1.25