This saves computation for M&M by precomputing and re-using quantitaties shared between iterations. It mostly saves $O(R^3)$ computations. This vignette shows results agree with the original version. Cannot use unit test due to numerical discrepency between chol of amardillo and R -- this has been shown problematic for some computations. I'll have to improve mashr for it.
muffled_chol = function(x, ...)
withCallingHandlers(chol(x, ...),
warning = function(w) {
if (grepl("the matrix is either rank-deficient or indefinite", w$message))
invokeRestart("muffleWarning")
} )
set.seed(1)
library(mashr)
simdata = simple_sims(500,5,1)
data = mash_set_data(simdata$Bhat, simdata$Shat, alpha = 0)
U.c = cov_canonical(data)
grid = mashr:::autoselect_grid(data,sqrt(2))
Ulist = mashr:::normalize_Ulist(U.c)
xUlist = expand_cov(Ulist,grid,TRUE)
llik_mat0 = mashr:::calc_lik_rcpp(t(data$Bhat),t(data$Shat),data$V,
matrix(0,0,0), simplify2array(xUlist),T,T)$data
svs = data$Shat[1,] * t(data$V * data$Shat[1,])
sigma_rooti = list()
for (i in 1:length(xUlist)) sigma_rooti[[i]] = t(backsolve(muffled_chol(svs + xUlist[[i]], pivot=T), diag(nrow(svs))))
llik_mat = mashr:::calc_lik_common_rcpp(t(data$Bhat),
simplify2array(sigma_rooti),
T)$data
head(llik_mat0)
head(llik_mat)
rows <- which(apply(llik_mat,2,function (x) any(is.infinite(x))))
if (length(rows) > 0)
warning(paste("Some mixture components result in non-finite likelihoods,",
"either\n","due to numerical underflow/overflow,",
"or due to invalid covariance matrices",
paste(rows,collapse=", "), "\n"))
loglik_null = llik_mat[,1]
lfactors = apply(llik_mat,1,max)
llik_mat = llik_mat - lfactors
mixture_posterior_weights = mashr:::compute_posterior_weights(1/ncol(llik_mat), exp(llik_mat))
post0 = mashr:::calc_post_rcpp(t(data$Bhat), t(data$Shat), matrix(0,0,0), matrix(0,0,0),
data$V,
matrix(0,0,0), matrix(0,0,0),
simplify2array(xUlist),
t(mixture_posterior_weights),
T, 4)
Vinv = solve(svs)
U0 = list()
for (i in 1:length(xUlist)) U0[[i]] = xUlist[[i]] %*% solve(Vinv %*% xUlist[[i]] + diag(nrow(xUlist[[i]])))
post = mashr:::calc_post_precision_rcpp(t(data$Bhat), t(data$Shat), matrix(0,0,0), matrix(0,0,0),
data$V,
matrix(0,0,0), matrix(0,0,0),
Vinv,
simplify2array(U0),
t(mixture_posterior_weights),
4)
head(post$post_mean)
head(post0$post_mean)
head(post$post_cov)
head(post0$post_cov)
Now test the relevant mvsusieR interface:
simulate_multivariate = function(n=100,p=100,r=2) {
set.seed(1)
res = mvsusieR::mvsusie_sim1(n,p,r,4,center_scale=TRUE)
res$L = 10
return(res)
}
attach(simulate_multivariate(r=2))
prior_var = V[1,1]
residual_var = as.numeric(var(y))
data = mvsusieR:::DenseData$new(X,y)
A = mvsusieR:::BayesianSimpleRegression$new(ncol(X), residual_var, prior_var)
A$fit(data, save_summary_stats = T)
null_weight = 0
mash_init = mvsusieR:::MashInitializer$new(list(V), 1, 1 - null_weight, null_weight)
residual_covar = cov(y)
mash_init$precompute_cov_matrices(data, residual_covar)
B = mvsusieR:::MashRegression$new(ncol(X), residual_covar, mash_init)
B$fit(data, save_summary_stats = T)