samplePoisMixture {telescope}

Telescoping sampling of a Bayesian finite Poisson mixture with a prior on the number of components K.

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) and Früwirth-Schnatter and Malsiner-Walli (2019), see the vignette for details and notation.

Usage

samplePoisMixture(
  y,
  S,
  mu,
  eta,
  a0,
  b0,
  h0,
  H0,
  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 rate values.

eta

A numeric vector; containing the initial cluster sizes.

a0

A numeric vector; hyperparameter of the prior on the rate \mu.

b0

A numeric vector; hyperparameter of the prior on the rate \mu.

h0

A numeric vector; hyperparameter of the prior on the rate \mu.

H0

A numeric vector; hyperparameter of the prior on the rate \mu.

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 rate \mu.

• "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 e0 or \alpha.

Examples

N <- 200
z <- sample(1:2, N, prob = c(0.5, 0.5), replace = TRUE)
y <- rpois(N, c(1, 6)[z])

Mmax <- 200
thin <- 1
burnin <- 100
M <- Mmax/thin

Kmax <- 50  
Kinit <- 10

G <- "MixDynamic"
priorOnAlpha <- priorOnAlpha_spec("gam_1_2")
priorOnK <- priorOnK_spec("BNB_143")

a0 <- 0.1 
h0 <- 0.5 
b0 <- a0/mean(y) 
H0 <- h0/b0

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)

result <- samplePoisMixture(
  y, S_0, mu_0, eta_0, 
  a0, b0, h0, H0,
  M, burnin, thin, Kmax, 
  G, priorOnK, priorOnAlpha)

K <- result$K
Kplus <- result$Kplus

plot(seq_along(K), K, type = "l", ylim = c(0, max(K)), 
     xlab = "iteration", main = "",
     ylab = expression("K" ~ "/" ~ K["+"]), col = 1)
lines(seq_along(Kplus), Kplus, col = 2)
legend("topright", legend = c("K", expression(K["+"])),
       col = 1:2, lty = 1, box.lwd = 0)