# R/VB_bms.R In mattelisi/bmsR: Bayesian model selection for group studies in R

#### Documented in VB_bms

```#' Bayesian model selection for group studies
#'
#' This function is translated from spm_BMS.m function, of SPM12, https://github.com/neurodebian/spm12/blob/master/spm_BMS.m
#' @param m N by K (N subjects by K models) matrix of log-model evidence.
#' @param n_samples Number of samples used in the simulation used for calculating exceedance probabilities.
#' @return a list with the posterior estimates of the Dirichlet parameters (alpha), the expected model frequencies (r), the exceedance probabilities (xp), the Bayesian Omnibus Risk (bor), and the protected exceedance probabilities (pxp).
#' @export

VB_bms <-function(m, n_samples=1e6){

max_val <- log(.Machine\$double.xmax)
Ni      <- dim(m)[1]  # number of subjects
Nk      <- dim(m)[2]  # number of models
c       <- 1
cc      <- 10e-4
# norm_vec <- function(x) sqrt(sum(x^2))

alpha0 <- rep(1, Nk)

log_u <- matrix(NA,dim(m)[1],dim(m)[2])
u <- matrix(NA,dim(m)[1],dim(m)[2])
g <- matrix(NA,dim(m)[1],dim(m)[2])
alpha <- alpha0
beta <- rep(NA,Nk)

# iterative VB estimation
while(c > cc){
for(i in 1:Ni){
for(k in 1:Nk){
# integrate out prior probabilities of models (in log space)
log_u[i,k] <- m[i,k] + digamma(alpha[k])- digamma(sum(alpha))
}
# prevent numerical problems for badly scaled posteriors
for(k in 1:Nk){
log_u[i,k] <- sign(log_u[i,k]) * min(max_val, abs(log_u[i,k]))
}

# exponentiate (to get back to non-log representation)
u[i,]  <- exp(log_u[i,])

# normalisation: sum across all models for i-th subject
u_i  <- sum(u[i,])
g[i,] <- u[i,]/u_i
}

# expected number of subjects whose data we believe to have been
# generated by model k
for(k in 1:Nk){
beta[k] <- sum(g[,k])
}

# update alpha
prev  <- alpha
for(k in 1:Nk){
alpha[k] <- alpha0[k] + beta[k]
}

# convergence?
c <- sqrt(sum((alpha - prev)^2))

}

# the posterior p(r|y)
exp_r <- alpha/sum(alpha)

# exceedance probabilities
xp <- dirichlet_exceedance(alpha,n_samples)

# Bayes Omnibus Risk and Preotected xp
F1 <- FE(m, alpha , g) # Evidence of alternative
F0 <- FE_null(m) # Evidence of null (equal model freqs)

# Implied by Eq 5 (see also p39) in Rigoux et al.