R/dirichlet_exceedance.R

Defines functions dirichlet_exceedance

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.