<jupyter_start><jupyter_text>Pre-computing various second-moment related quantitiesThis 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.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Now test the relevant `mmbr` interface:<jupyter_code>simulate_multivariate = function(n=100,p=100,r=2) {
  set.seed(1)
  res = mmbr::mmbr_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 = mmbr:::DenseData$new(X,y)
A = mmbr:::BayesianSimpleRegression$new(ncol(X), residual_var, prior_var)
A$fit(data, save_summary_stats = T)
null_weight = 0
mash_init = mmbr:::MashInitializer$new(list(V), 1, 1 - null_weight, null_weight)
residual_covar = cov(y)
mash_init$precompute_cov_matrices(data, residual_covar)
B = mmbr:::MashRegression$new(ncol(X), residual_covar, mash_init)
B$fit(data, save_summary_stats = T)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13739858675739336, 'fraction_numerical': 0.02015179272441769, 'mean_word_length': 3.1141011840688915, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '10140', 'license': 'MIT', 'n_tokens_mistral': 1418, 'n_tokens_neox': 1340, 'n_words': 239, 'path': 'inst/prototypes/precomputed_quantities_mash.ipynb', 'repo_name': 'zouyuxin/mmbr', 'size': '3821'}