gumbel_max: Simulate Cluster Labels from Unnormalised Log-Probabilities using the Gumbel-Max Trick

Description

Samples cluster labels for N observations from G groups efficiently using log-probabilities and the so-called Gumbel-Max trick, without requiring that the log-probabilities be normalised; thus redunant computation can be avoided. Computation takes place on the log scale for stability/underflow reasons (to ensure negligible probabilities won't cause computational difficulties); in any case, many functions for calculating multivariate normal densities already output on the log scale.

Usage

gumbel_max(probs, slice = FALSE)

Arguments

probs

An N x G matrix of unnormalised probabilities on the log scale, where N is he number of observations that require labels to be sampled and G is the number of active clusters s.t. sampled labels can take values in 1:G.

slice

A logical indicating whether or not the indicator correction for slice sampling has been applied to probs. Defaults to FALSE but is TRUE for the "IMIFA" and "IMFA" methods under mcmc_IMIFA. Details of this correction are given in Murphy et. al. (2017). When set to TRUE, this results in a speed-improvement when probs contains non-finite values (e.g. -Inf, corresponding to zero on the probability scale).

Value

A vector of N sampled cluster labels, with the largest label no greater than G.

References

Murphy, K., Gormley, I. C. and Viroli, C. (2017) Infinite Mixtures of Infinite Factor Analysers: Nonparametric Model-Based Clustering via Latent Gaussian Models, arXiv:1701.07010.

Yellot, J. I. Jr. (1977) The relationship between Luce's choice axiom, Thurstone's theory of comparative judgment, and the double exponential distribution, Journal of Mathematical Psychology, 15: 109-144.

Examples

Run this code

# NOT RUN {
# Create a 1-row matrix of weights
  G         <- 3
  weights   <- matrix(c(1, 2, 3), nrow=1, ncol=G)

# Call gumbel_max() repeatedly to obtain samples of the labels, zs
  iters     <- 10000
  zs        <- vapply(seq_len(iters), function(i)
               gumbel_max(probs=log(weights)), numeric(1L))

# Compare answer to the normalised weights
  tabulate(zs, nbins=G)/iters
  normalised <- as.numeric(weights/sum(weights))

# Simulate a matrix of dirichlet weights & the associated vector of N labels
  N       <- 400
  G       <- 8
  sizes   <- seq(from=85, to=15, by=-10)
  weights <- matrix(rDirichlet(N * G, alpha=1, nn=sizes), byrow=TRUE, nrow=N, ncol=G)
  zs      <- gumbel_max(probs=log(weights))
# }