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.
PYregression(y, x, mcmc = list(), prior = list(), output = list())a vector of observations, univariate dependent variable;
a vector of observations, univariate independent variable;
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.
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).
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 and the identity matrix);
a0 and b0 are the shape and scale parameters of the inverse gamma base measure on the scale component
(default are 2 and 1).
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 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 identity matrix); See details.
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 vector 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.
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.
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\) and 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).
Canale, A., Corradin, R., Nipoti, B. (2019), Importance conditional sampling for Bayesian nonparametric mixtures, arXiv preprint, arXiv:1906.08147
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
Kalli, M., Griffin, J. E., and Walker, S. G. (2011), Slice sampling mixture models. Statistics and Computing 21, 93-105.
Neal, R. M. (2000), Markov Chain Sampling Methods for Dirichlet Process Mixture Models, Journal of Computational and Graphical Statistics 9, 249-265.
# NOT RUN {
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)
# }
Run the code above in your browser using DataLab