In stephenslab/udr: Ultimate Deconvolution for Multivariate Normal Means

Description:

This file is used to produce results for each type of algorithm at the first iteration for benchmarking purposes.

Unconstrained subproblem

# Simulate data for unconstrained subproblem
seed <- 1
n <- 500
U <- diag(2)
V <- diag(2)
X <- simulate_one_component(seed, n, U, V)

# Initialization
p <- runif(n)
U.init <- sim_unconstrained(2)
U.init <- create_unconstrained_matrix_struct(X, U.init)

# Unconstrained updates: result for one iteration 
U_ted <- update_prior_covariance_unconstrained_ted(X, U.init, V, p, minval = 0)$mat
U_fa <- update_prior_covariance_unconstrained_fa(X, U.init$mat, p)
U_ed <- update_prior_covariance_ed_iid(X, U.init$mat,V, p)

save(U_ted, file = 'U_ted.rds')
save(U_fa, file = 'U_fa.rds')
save(U_ed, file = 'U_ed.rds')

Scaled subproblem

# Simulate data
seed <- 1
n <- 500
U0 <- diag(c(1,1,0))
V <- diag(3)
X <- simulate_one_component(seed, n, U0, V)

# Initialization
U <- create_scaled_matrix_struct(X, U0)
minval <- 0
p <- runif(n)
r <- sum(eigen(U$U0)$values > 1e-15)

scalar  <- update_prior_covariance_scaled_fa_iid(X, U$U0, p, U$s, r)
U_scalar <- update_prior_covariance_scaled_struct(U, scalar)
save(U_scalar, file = 'scalar.rds')

Rank1 subproblem

# Simulate data
seed <- 1
n <- 500
U.true <- matrix(c(1,0,0,0), ncol = 2, nrow = 2)
V <- diag(2)
V.g <- simplify2array(replicate(n, V, simplify=FALSE))
X <- simulate_one_component(seed, n, U.true, V)

# Intialization 
U <- create_rank1_matrix_struct(X, sim_rank1(2))
p <- runif(n)
minval <- 0

# Rank1 updates: result for one iteration 
rank1_ted <- update_prior_covariance_rank1_ted(X, U, V, p, minval)$vec
rank1_fa <- update_prior_covariance_rank1_fa_iid(X, U$vec, p)
save(rank1_ted, file = 'rank1_ted.rds')
save(rank1_fa, file = 'rank1_fa.rds')