# nealAlgorithm3: Conjugate Gibbs Sampler for a Partition In sams: Merge-Split Samplers for Conjugate Bayesian Nonparametric Models

## Description

Algorithm 3 from Neal (2000) to update the state of a partition based on the "Chinese Restaurant Process" (CRP) prior and a user-supplied log posterior predictive density function, with additional functionality for the two parameter CRP prior.

## Usage

 ```1 2 3 4 5 6 7``` ```nealAlgorithm3( partition, logPosteriorPredictiveDensity = function(i, subset) 0, mass = 1, discount = 0, nUpdates = 1L ) ```

## Arguments

 `partition` A numeric vector of cluster labels representing the current partition. `logPosteriorPredictiveDensity` A function taking an index i (as a numeric vector of length one) and a subset of integers subset, and returning the natural logarithm of p( y_i | y_subset ), i.e., that item's contribution to the log integrated likelihood given a subset of the other items. The default value "turns off" the likelihood, resulting in prior simulation (rather than posterior simulation). `mass` A specification of the mass (concentration) parameter in the CRP prior. Must be greater than the `-discount` argument. `discount` A numeric value on the interval [0,1) corresponding to the discount parameter in the two parameter CRP prior. Set to zero for the usual, one parameter CRP prior. `nUpdates` An integer giving the number of Gibbs scans before returning. This has the effect of thinning the Markov chain.

## Value

A numeric vector giving the updated partition encoded using cluster labels.

## References

Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. Journal of computational and graphical statistics, 9(2), 249-265.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28``` ```nealData <- c(-1.48, -1.40, -1.16, -1.08, -1.02, 0.14, 0.51, 0.53, 0.78) mkLogPosteriorPredictiveDensity <- function(data = nealData, sigma2 = 0.1^2, mu0 = 0, sigma02 = 1) { function(i, subset) { posteriorVariance <- 1 / ( 1/sigma02 + length(subset)/sigma2 ) posteriorMean <- posteriorVariance * ( mu0/sigma02 + sum(data[subset])/sigma2 ) posteriorPredictiveSD <- sqrt(posteriorVariance + sigma2) dnorm(data[i], posteriorMean, posteriorPredictiveSD, log=TRUE) } } logPostPredict <- mkLogPosteriorPredictiveDensity() nSamples <- 1000L partitions <- matrix(0, nrow = nSamples, ncol = length(nealData)) for (i in 2:nSamples) { partitions[i,] <- nealAlgorithm3(partitions[i-1,], logPostPredict, mass = 1.0, nUpdates = 1) } # convergence and mixing diagnostics nSubsets <- apply(partitions, 1, function(x) length(unique(x))) mean(nSubsets) sum(acf(nSubsets)\$acf) - 1 # Autocorrelation time entropy <- apply(partitions, 1, partitionEntropy) plot.ts(entropy) ```

sams documentation built on June 16, 2021, 1:06 a.m.