Gibbs Sampling

In Bayesian framework, unknown variables were sampled and updated from the conditional posterior distribution using Markov Chain Monte Carlo (MCMC) [14]. Considering the likelihood and priors in Formulae 3 and 4, the full joint posterior distribution can be written as follows:

Unobservables (c, μ, σ²) were repeatedly sampled and updated from their posteriors, conditional on all other variables. The Gibbs sampler was implemented as follows:

Initialization: Assign initial values for (μ_k, σ_k²) where k = 1 and c_i = 1, for i = 1 : n.

Update θ_i: The conditional posterior distribution of θ_i was

Update cluster indicators c_{i :} The conditional posterior probabilities for c_i were:

Pci=K+1|else∝α∬NyiθiσK+12NθiμK+1τi2NμK+1μ0σ02IGσK+12r1r2dμK+1dσK+12∝α∫NyiθiσK+12IGσK+12r1r2dσK+12∫NθiμK+1τi2NμK+1μ0σ02dμK+1=α2πr2r1Γr1Γr1+1212yi−θi2+r2r1+121τi2+σ02exp−θi−μ022τi2+σ02 where Γ(.) is the gamma function. Note that constant 1n−1+α was omitted in both probabilities and (μ_K + 1, σ_K + 1²) were unknown and needed to be integrated out to leave c_i as the only variable to be estimated from the Markov Chain. The Dirichlet Process was represented via the CRP [15]. Effects were assigned to either currently holding cluster(s) or a new cluster based on the above probabilities. If a new cluster was chosen, then the cluster size was increased, i.e. K + 1 → K. In case of n_− i,k = 0, the k^th cluster was eliminated and the cluster indicators were decreased by one, i.e. K → K − 1.

Resample and update (μ_k, σ_k²) suggested by Formulae 2 as per Neal [12] as follows:

where n_k is the number of effects associated with the k^th mixture component. The derivations of the fully conditional posterior distributions are detailed in the Appendix.

Repeat Steps 2 to 4.

Gibbs sampler was implemented with 100,000 iterations of the MCMC to update conditional posterior distributions. The first 80,000 samples were discarded as burn-in and the rest of the 20,000 samples were used to construct joint posterior distribution. The hyper-parameters in Algorithm 5 were set to be α = 0.05, r₁ = 1, r₂ = 0.01, μ₀ = 0, σ₀² = 0.01. Among hyperparameters, alpha was empirically set to 0.05 based on the simulation results. (Note: The larger the magnitude of alpha, the higher the probability of a large number of clusters.) Convergence was checked by inspection of negative log-likelihood plots. After the burn-in period, when the MCMC converges to the stationary distribution, sampled parameters were collected to form the posterior distribution. We employed posterior means for estimating the mean and variance μ^kσ^k2 and posterior modes for estimating ĉ_i, which was further used to infer π^k. The Bayesian confidence interval (BCI), which is the counterpart of the confidence interval in frequentist statistics, was defined as posterior probability that the parameter lies within the interval:

where α is the significance level. Instead of analytically estimating the confidence interval, the confidence interval for μ^kσ^k2 was numerically estimated from quartiles of posterior distribution.