sampleUniNormMixture: Telescoping sampling of a Bayesian finite univariate Gaussian...

In telescope: Bayesian Mixtures with an Unknown Number of Components

Telescoping sampling of a Bayesian finite univariate Gaussian mixture where a prior on the number of components K is specified.

Description

• The MCMC scheme is implemented as suggested in Frühwirth-Schnatter et al (2021).

• The priors on the model parameters are specified as in Frühwirth-Schnatter et al (2021), see the vignette for details and notation.

• The parametrizations of the gamma and inverse gamma distribution are used as in Frühwirth-Schnatter et al (2021), see also the vignette.

Usage

sampleUniNormMixture(
  y,
  S,
  mu,
  sigma2,
  eta,
  c0,
  g0,
  G0,
  C0_0,
  b0,
  B0,
  M,
  burnin,
  thin,
  Kmax,
  G = c("MixDynamic", "MixStatic"),
  priorOnK,
  priorOnWeights,
  verbose = FALSE
)

Arguments

y	A numeric matrix; containing the data.
S	A numeric matrix; containing the initial cluster assignments.
mu	A numeric matrix; containing the initial cluster-specific mean values.
sigma2	A numeric matrix; containing the initial cluster-specific variance values.
eta	A numeric vector; containing the initial cluster sizes.
c0	A numeric vector; hyperparameter of the prior on \sigma^2_k.
g0	A numeric vector; hyperparameter of the prior on \sigma^2_k
G0	A numeric vector; hyperparameter of the prior on \sigma^2_k
C0_0	A numeric vector; initial value of hyperparameter C_0.
b0	A numeric vector; hyperparameter of the prior on \mu_k.
B0	A numeric vector; hyperparameter of the prior on \mu_k.
M	A numeric scalar; specifying the number of recorded iterations.
burnin	A numeric scalar; specifying the number of burn-in iterations.
thin	A numeric scalar; specifying the thinning used for the iterations.
Kmax	A numeric scalar; the maximum number of components.
G	A character string; either "MixDynamic" or "MixStatic".
priorOnK	A named list; providing the prior on the number of components K, see priorOnK_spec().
priorOnWeights	A named list; providing the prior on the mixture weights.
verbose	A logical; indicating if some intermediate clustering results should be printed.

Value

A named list containing:

• "Mu": sampled component means.

• "Sigma2": sampled component component variances.

• "Eta": sampled weights.

• "S": sampled assignments.

• "Nk": number of observations assigned to the different components, for each iteration.

• "K": sampled number of components.

• "Kplus": number of filled, i.e., non-empty components, for each iteration.

• "e0": sampled Dirichlet parameter of the prior on the weights (if e_0 is random).

• "alpha": sampled Dirichlet parameter of the prior on the weights (if \alpha is random).

• "acc": logical vector indicating acceptance in the Metropolis-Hastings step when sampling either e_0 or \alpha.

Examples

if (requireNamespace("mclust", quietly = TRUE)) {
    data("acidity", package = "mclust")
    y <- acidity
    
    N <- length(y)
    r <- 1
    
    M <- 200
    thin <- 1
    burnin <- 100
    Kmax <- 50  
    Kinit <- 10
    
    G <- "MixStatic" 
    priorOnE0 <- priorOnE0_spec("e0const", 0.01)
    priorOnK <- priorOnK_spec("Pois_1", 50)
    
    R <- diff(range(y))
    c0 <- 2 + (r-1)/2
    C0 <- diag(c(0.02*(R^2)), nrow = r)
    g0 <- 0.2 + (r-1) / 2
    G0 <- diag(10/(R^2), nrow = r)
    B0 <- diag((R^2), nrow = r)
    b0 <- as.matrix((max(y) + min(y))/2, ncol = 1)  
    
    cl_y <- kmeans(y, centers = Kinit, nstart = 100)
    S_0 <- cl_y$cluster
    mu_0 <- t(cl_y$centers)
    eta_0 <- rep(1/Kinit, Kinit)
    sigma2_0 <- array(0, dim = c(1, 1, Kinit))
    sigma2_0[1, 1, ] <- 0.5 * C0

    result <- sampleUniNormMixture(
        y, S_0, mu_0, sigma2_0, eta_0,
        c0, g0, G0, C0, b0, B0,
        M, burnin, thin, Kmax,
        G, priorOnK, priorOnE0)
    
    K <- result$K
    Kplus <- result$Kplus  
    
    plot(K, type = "l", ylim = c(0, max(K)),
         xlab = "iteration", main = "",
         ylab = expression("K" ~ "/" ~ K["+"]), col = 1)
    lines(Kplus, col = 2)
    legend("topright", legend = c("K", expression(K["+"])),
           col = 1:2, lty = 1, box.lwd = 0)
}

telescope documentation built on April 4, 2025, 2:40 a.m.