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

#### Documented in dirichlet_exceedance

```#' Compute exceedance probabilities
#'
#' This function compute exceedance probabilities by simulating the multivariate Dirichlet observations using the same approach implemented in SPM 12.
#' @param alpha Dirichlet parameters (usually the posterior estimate).
#' @param Nsamp number of samples drawn.
#' @return Exceedance probabilities.
#' @export

dirichlet_exceedance <- function(alpha, Nsamp=1e6){
Nk <- length(alpha)
ind_max <- function(x) ifelse(x==max(x),1,0)
xp <- rep(0,Nk)

# Sample from univariate gamma densities then normalise
# (see Dirichlet entry in Wikipedia or Ferguson (1973) Ann. Stat. 1,
#   % 209-230)
r = mlisi::repmat(as.matrix(0), Nsamp, Nk)

for(k in 1:Nk){
r[,k] <- rgamma(Nsamp, alpha[k], scale = 1)
}
sr <- apply(r,1,sum)

for(k in 1:Nk){
r[,k] <- r[,k]/sr
}

# Exceedance probabilities:
# For any given model k1, compute the probability that it is more
# likely than any other model k2~=k1
j <- apply(r, 1, which.max)
xp <- xp + unname(tapply(j, factor(j),length))
xp <-xp/Nsamp
}
```
mattelisi/bmsR documentation built on Jan. 29, 2019, 6:20 p.m.