PYregression {BNPmix}

MCMC for Pitman-Yor mixture of Gaussian regressions

Description

The PYregression function generates a posterior sample for mixtures of linear regression models inspired by the ANOVA-DDP model introduced in De Iorio et al. (2004). See details below for model specification.

Usage

PYregression(y, x, mcmc = list(), prior = list(), output = list())

Arguments

y

a vector of observations, univariate dependent variable;

x

a matrix of observations, multivariate independent variable;

mcmc

a list of MCMC arguments:

• niter (mandatory), number of iterations.

• nburn (mandatory), number of iterations to discard as burn-in.

• method, the MCMC sampling method to be used. Options are 'ICS', 'MAR' and 'SLI' (default is 'ICS'). See details.

• model the type of model to be fitted (default is 'LS'). See details.

• nupd, argument controlling the number of iterations to be displayed on screen: the function reports on standard output every time nupd new iterations have been carried out (default is niter/10).

• print_message, control option. If equal to TRUE, the status is printed to standard output every nupd iterations (default is TRUE).

• m_imp, number of generated values for the importance sampling step of method = 'ICS' (default is 10). See details.

• slice_type, when method = 'SLI' it specifies the type of slice sampler. Options are 'DEP' for dependent slice-efficient, and 'INDEP' for independent slice-efficient (default is 'DEP'). See details.

• m_marginal, number of generated values for the augmentation step needed, if method = 'MAR', to implement Algorithm 8 of Neal, 2000. (Default is 100). See details.

• hyper, if equal to TRUE, hyperprior distributions on the base measure's parameters are added, as specified in prior and explained in details (default is TRUE).

prior

a list giving the prior information. The list includes strength and discount, the strength and discount parameters of the Pitman-Yor process (default are strength = 1 and discount = 0, the latter leading to the Dirichlet process). The remaining parameters specify the base measure: m0 and S0 are the mean and covariance of normal base measure on the regression coefficients (default are a vector of zeroes, except for the first element equal to mean(y), and a diagonal matrix with each element equal to 100); a0 and b0 are the shape and scale parameters of the inverse gamma base measure on the scale component (default are 2 and var(y)). If hyper = TRUE, optional hyperpriors on the base measure's parameters are added: specifically, m1 and k1 are the mean parameter and scale factor defining the normal hyperprior on m0 (default are a vector of zeroes, except for the first element equal to the sample mean of the dependent observed variable, and 1); tau1 and zeta1 are the shape and rate parameters of the gamma hyperprior on b0 (default is 1 for both); n1 and S1 are the parameters (degrees of freedom and scale) of the Wishart prior for S0 (default 4 and a diagonal matrix with each element equal to 100); See details.

output

list of posterior summaries:

• grid_y, a vector of points where to evaluate the estimated posterior mean density of y, conditionally on each value of x in grid_x;

• grid_x, a matrix of points where to evaluate the realization of the posterior conditional densities of y given x;

• out_type, if out_type = "FULL", the function returns the estimated partitions and the realizations of the posterior density for each iteration; If out_type = "MEAN", return the estimated partitions and the mean of the densities sampled at each iteration; If out_type = "CLUST", return the estimated partitions. Default out_type = "FULL";

• out_param, if equal to TRUE, the function returns the draws of the kernel's parameters for each MCMC iteration, default is FALSE. See value for details.

Details

This function fits a Pitman-Yor process mixture of Gaussian linear regression models, i.e

\tilde f(y) = \int \phi(y; x^T \beta, \sigma^2) \tilde p (d \beta, d \sigma^2)

where x is a bivariate vector containing the dependent variable in x and a value of 1 for the intercept term. The mixing measure \tilde p has a Pitman-Yor process prior with strength \vartheta, discount parameter \alpha. The location model assume a base measures P_0 specified as

P_0(d \beta) = N(d \beta; m_0, S_0) .

while the location-scale model assume a base measures P_0 specified as

P_0(d \beta, d \sigma^2) = N(d \beta; m_0, S_0) \times IGa(d \sigma^2; a_0, b_0).

Optional hyperpriors complete the model specification:

m_0 \sim N(m_1, S_0 / k_1 ),\quad S_0 \sim IW(\nu_1, S_1),\quad b_0 \sim G(\tau_1, \zeta_1).

Posterior simulation methods

This generic function implements three types of MCMC algorithms for posterior simulation. The default method is the importance conditional sampler 'ICS' (Canale et al. 2019). Other options are the marginal sampler 'MAR' (algorithm 8 of Neal, 2000) and the slice sampler 'SLI' (Kalli et al. 2011). The importance conditional sampler performs an importance sampling step when updating the values of individual parameters \theta, which requires to sample m_imp values from a suitable proposal. Large values of m_imp are known to improve the mixing of the posterior distribution at the cost of increased running time (Canale et al. 2019). When updateing the individual parameter \theta, Algorithm 8 of Neal, 2000, requires to sample m_marginal values from the base measure. m_marginal can be chosen arbitrarily. Two options are available for the slice sampler, namely the dependent slice-efficient sampler (slice_type = 'DEP'), which is set as default, and the independent slice-efficient sampler (slice_type = 'INDEP') (Kalli et al. 2011). See Corradin et al. (to appear) for more details.

Value

A BNPdens class object containing the estimated density and the cluster allocations for each iterations. The output contains also the data and the grids. If out_param = TRUE the output contains also the kernel specific parameters for each iteration. If mcmc_dens = TRUE, the function returns also a realization from the posterior density for each iteration. If mean_dens = TRUE, the output contains just the mean of the densities sampled at each iteration. The output retuns also the number of iterations, the number of burn-in iterations, the computational time and the type of model.

References

Canale, A., Corradin, R., Nipoti, B. (2019), Importance conditional sampling for Bayesian nonparametric mixtures, arXiv preprint, arXiv:1906.08147

Corradin, R., Canale, A., Nipoti, B. (2021), BNPmix: An R Package for Bayesian Nonparametric Modeling via Pitman-Yor Mixtures, Journal of Statistical Software, doi:10.18637/jss.v100.i15

De Iorio, M., Mueller, P., Rosner, G.L., and MacEachern, S. (2004), An ANOVA Model for Dependent Random Measures, Journal of the American Statistical Association 99, 205-215, doi:10.1198/016214504000000205

Kalli, M., Griffin, J. E., and Walker, S. G. (2011), Slice sampling mixture models. Statistics and Computing 21, 93-105, doi:10.1007/s11222-009-9150-y

Neal, R. M. (2000), Markov Chain Sampling Methods for Dirichlet Process Mixture Models, Journal of Computational and Graphical Statistics 9, 249-265, doi:10.2307/1390653

Examples

x_toy <- c(rnorm(100, 3, 1), rnorm(100, 3, 1))
y_toy <- c(x_toy[1:100] * 2 + 1, x_toy[101:200] * 6 + 1) + rnorm(200, 0, 1)
grid_x <- c(0, 1, 2, 3, 4, 5)
grid_y <- seq(0, 35, length.out = 50)
est_model <- PYregression(y = y_toy, x = x_toy,
mcmc = list(niter = 200, nburn = 100),
output = list(grid_x = grid_x, grid_y = grid_y))
summary(est_model)
plot(est_model)