In stephenslab/mvebnm: Ultimate Deconvolution for Multivariate Normal Means

knitr::opts_chunk$set(comment = "#",collapse = TRUE,results = "hold",
                      fig.align = "center")

Here we illustrate how to incorporate covariance regularization to unconstrained update. This is recommended when sample size is small and the number of condition is large.

library(udr)

Simulate data

set.seed(3)
n = 20
R = 5
K = 2
Ut = c()
for (k in 1:K){
  Ut[[k]] <- udr:::sim_unconstrained(R)
}
Vt = diag(R)
wt = rep(1/K, K)
# Simulate data 
X <- simulate_ud_data(n,wt,Ut,Vt)
X = as.matrix(X)

Initialization

set.seed(999)
U.init = c()
for (i in 1:K)
  U.init[[i]] = udr:::sim_unconstrained(R)

tol = 1e-2
tol.lik = 1e-2
lambda = 1

f0 <- ud_init(X, V = Vt, n_rank1 = 0, U_scaled = NULL, U_unconstrained = U.init)

Run udr with covariance regularization

In control list, specify (1) $\lambda$: the parameter that controls regularization strength; (2) penalty.type: either "iw" or "nu", where "iw" means inverse Wishart penalty and "nu" means nuclear norm penalty. In total, there are three algorithms that can be used:

1. ted.nu: TED algorithm with nuclear norm penalty;

2. ted.iw: TED algorithm with inverse Wishart penalty;

3. ed.iw: ED algorithm with inverse Wishart penalty.

1. Run ted.nu

ted.nu = ud_fit(f0, control = list(unconstrained.update = "ted", scaled.update = "fa", resid.update = 'none', lambda =lambda, penalty.type = "nu", maxiter=500, tol = tol, tol.lik = tol.lik), verbose=TRUE)

2. Run ted.iw

ted.iw = ud_fit(f0, control = list(unconstrained.update = "ted", scaled.update = "fa", resid.update = 'none', lambda =lambda, penalty.type = "iw", maxiter=500, tol = tol, tol.lik = tol.lik), verbose=TRUE)

3. Run ed.iw

ed.iw = ud_fit(f0, control = list(unconstrained.update = "ed", scaled.update = "fa", resid.update = 'none', lambda =lambda, penalty.type = "iw", maxiter=500, tol = tol, tol.lik = tol.lik), verbose=TRUE)

Plot results:

The penalized log-likelihood is guaranteed to increase every iteration.

plot(ted.nu$progress$loglik.pen, ylim = c(-200, -188), xlim = c(0, 20), xlab = "iteration", ylab = "Penalized log-likelihood")
points(ted.iw$progress$loglik.pen, col = 2)
points(ed.iw$progress$loglik.pen, col = 3)
legend("bottomright", legend = c("ted.nu", "ted.iw", "ed.iw"), col = c(1:3), pch = 1)

stephenslab/mvebnm documentation built on June 4, 2024, 6:37 a.m.