A Conditional Density Estimation Partition Model Using Logistic Gaussian Processes

TL;DR

This paper introduces a novel conditional density estimation method using logistic Gaussian processes within a partition model framework, employing Voronoi tessellations and reversible jump MCMC for efficient inference.

Contribution

It presents a new approach combining logistic Gaussian processes with Voronoi tessellations and Laplace approximations, enabling effective conditional density estimation with proven consistency.

Findings

Successfully estimates partition structure and conditional distribution in simulations

Demonstrates effectiveness in real data applications

Provides a computationally feasible inference algorithm

Abstract

Conditional density estimation (density regression) estimates the distribution of a response variable y conditional on covariates x. Utilizing a partition model framework, a conditional density estimation method is proposed using logistic Gaussian processes. The partition is created using a Voronoi tessellation and is learned from the data using a reversible jump Markov chain Monte Carlo algorithm. The Markov chain Monte Carlo algorithm is made possible through a Laplace approximation on the latent variables of the logistic Gaussian process model. This approximation marginalizes the parameters in each partition element, allowing an efficient search of the posterior distribution of the tessellation. The method has desirable consistency properties. In simulation and applications, the model successfully estimates the partition structure and conditional distribution of y.

Equations

$$p(M)$$

$$p(\mathbf{c}\mid M)$$

$$p(\mathbf{w})$$

$$p(y_i)=\frac{\exp (f_i(y_i))}{{\int }_{{\mathcal{V}}_i}\exp (f_i(s))ds}$$

$$f_i(z_{ij})$$

$$g(z_{ij})$$

$${\mu }_{ij}$$

$$f_i(z_{ij})\sim GP(\mathbf{h}(z_{ij}{)}^T\mathbf{b},{\kappa }_{{\mathbf{\theta }}_i}(z_{ij},z_{i{j}^{{}^{\prime }}})+\mathbf{h}(z_{ij}{)}^{T}B\mathbf{h}(z_{i{j}^{{}^{\prime }}})).$$

$$p({\mathbf{f}}_i\mid Z_i,{\mathbf{\theta }}_i)=N({\mathbf{f}}_i\mid H_i\mathbf{b},K_i+H_{i}B{H}_i^T)$$

$$\log [p({\mathbf{y}}_i\mid {\mathbf{f}}_i,Z_i)]={\mathbf{y}}_i^{\star T}{\mathbf{f}}_i-n_i\log \left(\sum _{j=1}^r\exp ({\mathbf{f}}_{ij})\right)$$

$$p({\mathbf{f}}_i\mid {\mathbf{y}}_i,Z_i,{\mathbf{\theta }}_i)=\frac{p({\mathbf{f}}_i\mid Z_i,{\mathbf{\theta }}_i)p({\mathbf{y}}_i\mid {\mathbf{f}}_i,Z_i)}{{\int }_{{\mathcal{V}}_i}p(\mathbf{f}\mid Z_i,{\mathbf{\theta }}_i)p({\mathbf{y}}_i\mid \mathbf{f},Z_i)d\mathbf{f}}$$

$$p({\mathbf{y}}_i\mid Z_i,{\mathbf{\theta }}_i)={\int }_{{\mathcal{V}}_i}p({\mathbf{f}}_i\mid Z_i,{\mathbf{\theta }}_i)p({\mathbf{y}}_i\mid {\mathbf{f}}_i,Z_i)d{\mathbf{f}}_i.$$

$$p(T\mid \mathbf{y},\mathbf{\theta },[Z_i{]}_{i=1}^M)\propto p(T)\prod _{i=1}^{M}p({\mathbf{y}}_i\mid {\mathbf{\theta }}_i,Z_i,T).$$

$$p({\mathbf{f}}_i\mid {\mathbf{y}}_i,Z_i,{\mathbf{\theta }}_i)\approx N({\mathbf{f}}_i\mid {\hat{\mathbf{f}}}_i,{\Sigma }_i)$$

$${\kappa }_{{\mathbf{\theta }}_i}(z,z^{\prime })={\sigma }_i^2\exp \left(-\frac{1}{2l_i^2}(z-z^{\prime }{)}^2\right)$$

$$\alpha =\min \left(1,\frac{q(\mathbf{w}\mid {\mathbf{w}}^{(p)}){\prod }_{i=1}^{M_p}p({\mathbf{y}}_i^{\star }\mid T_p,Z,{\mathbf{\theta }}_i)}{q({\mathbf{w}}^{(p)}\mid \mathbf{w}){\prod }_{i=1}^{M}p({\mathbf{y}}_i^{\star }\mid T_c,Z,{\mathbf{\theta }}_i)}\right).$$

$$g_{\mathbf{x}}(y)$$

$$f(y\mid \mathbf{x})$$

$$N_1=\{\Vert {\eta }_j(y)-{\eta }_j^{\ast }(y){\Vert }_{\infty }<{ϵ}_2\text{ in }\mathbf{x}\in U_j^{\ast }\cap V_j^{\ast }\text{ and }\Vert {\eta }_j(y){\Vert }_{\infty }<k_0\text{ for }\mathbf{x}\in U_j^{\ast }\Delta V_j^{\ast }\},$$

$$N_2=\{\Vert {\mu }_j(y)-{\mu }_j^{\ast }(y){\Vert }_{\infty }<{ϵ}_2\text{ in }\mathbf{x}\in U_j^{\ast }\cap V_j^{\ast }\text{ and }\Vert {\mu }_j(y){\Vert }_{\infty }<k_0\text{ for }\mathbf{x}\in U_j^{\ast }\Delta V_j^{\ast }\}.$$

$$y\sim N(5,.5^2),x_1\sim N(0,1),x_2\sim N(0,5^2)$$

$$f(x)=\{\begin{array}{ll}2.5^2& x<2.5\\ x^2& x\ge 2.5\end{array}$$

$$x_1,x_2$$

$$y$$

$$y\mid x_1,x_2\sim \{\begin{array}{ll}N(3,.5^2)& x_1<5\\ -(x_1-.5{)}^2+Z_{x_2}& x_1\ge 5\end{array}$$

$$log(N(ϵ,{\mathbf{K}}_{j,\beta }^n,\Vert {\Vert }_{\infty }))=log(N(ϵ,{\mathbf{K}}_j^n,\Vert {\Vert }_{\infty }))+log(N(ϵ,{\mathbf{K}}_b^n,\Vert {\Vert }_{\infty }))=o(n).$$

$$log(\Pi (({\mathbf{K}}_{j,\beta }^n{)}^c)$$

$$log(N(ϵ,{\mathbf{K}}_{j,\beta }^n,\Vert {\Vert }_1))$$

$$log(\Pi ({\mathcal{K}}^c))=O(-n+logn)=O(-n)$$

$$log(N(ϵ,\mathcal{K},\Vert {\Vert }_1))=o(n+logn)=o(n),$$

$$\Pi (U_{{ϵ}^{\prime }}^c\mid \cdot )\le \frac{{\int }_{U_{{ϵ}^{\prime }}^c}{e}^{-log(f^{\ast }(\mathbf{y}))+log(f(\mathbf{y}))}d\Pi ()}{(1/M_{max}){\int }_{N_{1,2}^{{ϵ}_2}}{e}^{-log(f^{\ast }(\mathbf{y}))+log(f(\mathbf{y}))-Mlog(n)}d\Pi ()}=\frac{{\mathbf{N}}_1^n}{{\mathbf{D}}_1^n}.$$

$$\frac{{\mathbf{N}}_1^n}{{\mathbf{D}}_1^n}=e^{2n{\alpha }_1{ϵ}_3}\frac{{{\mathbf{N}}_1^n}_{\mathcal{K}}+{{\mathbf{N}}_1^n}_{{\mathcal{K}}^c}}{e^{2n{\alpha }_1{ϵ}_3}{\mathbf{D}}_1^n}$$

$$\Pi (U_{{ϵ}^{\prime }}^c\mid \cdot )={\Psi }_n\Pi (U_{{ϵ}^{\prime }}^c\mid \cdot )+(1-{\Psi }_n)\Pi (U_{{ϵ}^{\prime }}^c\mid \cdot )\le {\Psi }_n\Pi (U_{{ϵ}^{\prime }}^c\mid \cdot )+(1-{\Psi }_n)\frac{{\mathbf{N}}_1^n}{{\mathbf{D}}_1^n}.$$

$$E(\Pi (U_{{ϵ}^{\prime }}^c\mid \cdot ))\le e^{-c_{1}n{{ϵ}^{\prime }}^2/8}+e^{-c_{1}n{{ϵ}^{\prime }}^2/16}.$$

Full text

1 Introduction

Conditional density estimation (sometimes referred to as density regression) is a method used to estimate the conditional distribution of a response variable, $y$, which depends on a vector of covariates, $\mathbf{x}$. Many common regression methods are special cases of conditional density estimation. For instance, the Gaussian regression model assumes that the mean of $y$ changes with $\mathbf{x}$ with the variance of $y$ being constant over the covariate space. Conditional density estimation is particularly useful when a parametric form linking the covariate space $\mathbf{x}$ with $y$ is unknown or violates the assumptions of existing parametric methods. In its most flexible forms, conditional density estimation can be viewed as a general nonparametric regression method.

There are a number of existing frequentist approaches to perform density regression including kernel methods (Fan et al., 1996; Fu et al., 2011), spline methods (Kooperberg and Stone, 1991; Stone et al., 1997), and mixtures of experts (Jacobs et al., 1991). There are also several Bayesian approaches to conditional density estimation. One of the popular approaches is to use mixture models for the conditional distribution of $p(y\mid\mathbf{x})$ and allow the mixing weights as well as the parameters to depend on the covariates (Chung and Dunson, 2009; Dunson and Park, 2008; Dunson et al., 2007; Griffin and Steel, 2006). An alternative approach is to apply the logistic Gaussian process model in the conditional density estimation setting (Tokdar et al., 2010). Latent variable models have been utilized by Kundu and Dunson (2011) and Bhattacharya and Dunson (2010). A multivariate spline based method (Shen et al., 2016) and an optional Polya tree based method (Ma and Wong, 2011) have been recently proposed. Petralia et al. (2013) perform density regression using a convex combination of dictionary densities using a fixed tree decomposition which scales to accommodate a large number of features.

The purpose of this paper is to provide a novel partition model (Holmes et al., 2005; Denison et al., 2002) framework to perform density regression using logistic Gaussian processes. This proposed method is intended to describe the distribution of $y$ in a region of $\mathbf{x}$. This method is similar in spirit to the Bayesian classification and regression tree (CART) model (Chipman et al., 1998; Denison et al., 1998) as a decision tool to understand how and where the density of $y$ changes in different regions of $\mathbf{x}$. In fact, we adaptively partition the covariate space and use a logistic Gaussian process model within each region. This method can also be viewed as a non-parametric changepoint analysis in which the changepoints occur in a multi-dimensional covariate space and the distribution of $y$ is not restricted to any specific parametric form. Hence, we can derive a decision rule based on the partitions which makes it practically useful. For example, a business may be interested in understanding how years of experience and a persuasiveness score influence revenue from individual salespersons (where variability and skewness of revenue may change throughout the covariate space). The model can identify where major changes in the shape, spread, and/or center of the distribution of revenue occurs among their sales force. Another example is that of power output from windmills where assuming a fixed error structure of power output, $y$, over the covariate space (wind speed, direction, density, etc.) is not appropriate, nor is the relationship between $y$ and the covariates believed to be linear. While existing conditional density estimation methodologies allow the estimation of $y$ at a specific point, this methodology allows one to identify where the important changepoints exist, and estimate the density of $y$ in each region of $\mathbf{x}$. In addition, variable selection is automatically incorporated as irrelevant covariates will not influence the partitioning in this framework.

Partition models provide an appealing framework to determine the breakpoints in $\mathbf{x}$ where the density of $y$ changes. Specifically, we can utilize existing density estimation methods in constructing a density regression methodology. In the simplest cases, partition models fit independent/separate models to each piece of the partition. Ideally these partitions are determined using a data-driven approach. In our current endeavor, we aim to find a data-driven partition of the covariate space $\mathbf{x}$ and fit independent density estimates for $y$ in each partition. This type of partition methodology has been used in spatial applications by Kim et al. (2005) and in disease mapping applications by Denison and Holmes (2001).

We use a logistic Gaussian process model (Lenk, 1988, 1991) for conditional density estimation within each region of the partition. Tokdar et al. (2010) used similar models, but our method is completely different than their approach. They used the logistic Gaussian process to model the joint distribution of the response $y$ and the covariates $\mathbf{x}$ and utilized a subspace projection method to reduce the dimension of the covariates. On the other hand, we are fitting univariate logistic Gaussian processes within each region of $\mathbf{x}$, hence avoiding the curse of dimensionality. Furthermore, in joint modeling approaches, $y$ and $\mathbf{x}$ are not identified as the response and covariates. Our experience is that the distribution of the response $y$ is highly influenced by the distribution of $\mathbf{x}$ (especially when the dimension of $\mathbf{x}$ is high). In addition, the proposed method is highly interpretable as it explores the relationship between the covariates and the density function in a nonlinear way which may not be available in projection based methods. To the authors’ knowledge, there is no publicly available code to perform Bayesian conditional density estimation. Therefore, part of our contribution is providing the code to perform Bayesian density regression in the partition model framework. This code is publicly available at https://github.com/gitrichhub/bayes-cde.

Our proposed method is also different than the existing partition model approaches. The effectiveness of the partition model algorithms depends on explicit marginalization of the model parameters so that an efficient reversible jump Markov chain Monte Carlo (MCMC) algorithm (Green, 1995) can be developed over the number and location of the partitions. Therefore, all the existing partition models are based on the conjugate structure of the likelihood and priors. On the contrary, the logistic Gaussian model does not have a conjugate structure, so explicit marginalization is not possible. Hence, we utilize a Laplace approximation of the logistic Gaussian models (Riihimäki et al., 2014) to obtain the marginal likelihood in order to develop an efficient reversible jump MCMC algorithm. Furthermore, to the best of our knowledge, none of the existing papers investigated theoretical properties of the partition models. Ours is thus the first paper that considers the posterior consistency in estimating conditional distributions in the partition model framework. Indeed, there are a few papers which have considered theoretical properties of density regression models using other modeling frameworks (Tokdar and Ghosh, 2007; Pati et al., 2013; Norets and Pelenis, 2012; Bhattacharya and Dunson, 2010).

In Section 2 we present the conditional density estimation model in a partition framework and provide a reversible jump MCMC algorithm for estimation. Section 3 discusses some results on consistency. Section 4 applies the method to both simulated and real datasets and Section 5 concludes.

2 Bayesian Hierarchical Conditional Density Estimation Partition Model

2.1 Modeling the partition structures using a Voronoi tessellation

The partition model divides (partitions) the $p$-dimensional covariate space $\mathcal{D}$ into $M$ distinct pieces where $y$ is assumed to independently follow a different density $p_{i}(\cdot)$ within each partition. The partitioning of the covariate space is done through a Voronoi tessellation. The tessellation is defined by $M$ centers $\mathbf{c}_{1},\ldots,\mathbf{c}_{M}$ that divide the covariate space into $M$ disjoint regions $R_{1},\ldots,R_{M}$ where $R_{i}$ consists of all the observed $\mathbf{x}$ that are closest to center $\mathbf{c}_{i}$. Formally, $R_{i}=\{\mathbf{x}\in\mathcal{D}:||\mathbf{x}-\mathbf{c}_{i}||<||\mathbf{x}-\mathbf{c}_{j}||\ \forall\ i\neq j\}$. Here, $||\mathbf{x}||=||(x_{1},\ldots,x_{p})||=\sum_{i=1}^{p}w_{i}x_{i}^{2}$ where $\mathbf{w}$ is a normalized weighting vector ($\sum w_{k}=1$) which places different weights on each of the covariates (Holmes et al., 2005). The weighting provides additional flexibility in the tessellation and also performs variable selection (which will be demonstrated in the examples). Figure 1 shows an example of a Voronoi tessellation in two dimensions with $\mathbf{w}=(.5,.5)$.

For simplicity, we assume that the possible centers of the tessellation are restricted to the observed data points $\mathbf{x}$. Next, we assign prior distributions for both the number of centers as well as the center locations. An intuitive way to express the prior is $p(\mathbf{c},M,\mathbf{w})=p(\mathbf{c}\mid M)p(M)p(\mathbf{w})$ where

[TABLE]

where DU$(x\mid 1,\ldots,n)$ means discrete uniform on $1,\ldots,n$ and $M_{\text{max}}$ is the maximum number of allowable centers (a hyperparameter chosen by the user). The prior on $p(\mathbf{c}\mid M)$ gives equal weight to all possible combinations of $M$ centers with possible center locations corresponding to the $n$ observed values of the covariates $\mathbf{x}$. The vector $\mathbf{w}$ has the Dirichlet prior which is uniform on the simplex.

2.2 The data generating model

In this section, we develop the likelihood function by fitting a logistic Gaussian process model within each element of the partition. We assume that within the $i$th region the data follow a logistic-Gaussian model. Let $T=\{M,\mathbf{c},\mathbf{w}\}$ denote the tessellation parameters and $\bm{\theta}=(\bm{\theta}_{1},\ldots,\bm{\theta}_{M})$ denote the set of all other (smoothing) parameters in each partition region. Also let $\mathbf{y}_{i}=(y_{i1},\cdots,y_{in_{i}})$ and $\mathbf{x}_{i}$ denote the observed response variable and covariates found in the $i$th region of the partition, $i=1,\ldots,M$. Then the logistic Gaussian process model within the $i$th region will be

[TABLE]

where $f_{i}(y_{i})$ is an unconstrained latent function and $y_{i}$ lies in a finite region $\mathcal{V}_{i}$ of $\mathbb{R}$. A Gaussian process (GP) prior is placed on $f_{i}$, therefore smoothing can be controlled via the covariance structure.

We first discretize the region $\mathcal{V}_{i}$ into a regular grid with $r$ subregions centered at $Z_{i}=[z_{i1},\ldots,z_{ir}]^{T}$. We then specify the prior on $f_{i}$ by assuming

[TABLE]

where $\kappa_{\bm{\theta}_{i}}(z_{ij},z_{ij^{{}^{\prime}}})$ is the covariance function with inputs $z_{ij}$ and $z_{ij^{{}^{\prime}}}$ which depends on the smoothing hyperparameter $\bm{\theta}_{i}$. For simplicity in the following exposition, we assume $\bm{\theta}_{i}$ is fixed and known, and will discuss its selection later in this section. By placing a Gaussian prior on $\bm{\beta}_{i}\sim N(\mathbf{b},B)$, $\bm{\beta}_{i}$ can be integrated out to yield

[TABLE]

Consequently, given a discretization $Z_{i}$ and covariance parameters $\bm{\theta}_{i}$, the prior on $f_{i}$ can be expressed as a Gaussian distribution over the latent function values $\mathbf{f}_{i}$

[TABLE]

where $\mathbf{f}_{i}=[f_{i}(z_{ij})]_{j=1}^{r}$ is a column vector of $r$ latent values associated with each subregion, $K_{i}=[\kappa_{\bm{\theta}}(z_{ij},z_{ik})]_{j,k=1}^{r}$ and $H_{i}=[\mathbf{h}(z_{ij})^{T}]_{j=1}^{r}$. Following Riihimäki et al. (2014), we choose $\mathbf{h}(z_{ij})=[x,x^{2}]^{T}$.

Let $\mathbf{y}_{i}=(y_{i1},\cdots,y_{in_{i}})$ be the vector of $n_{i}$ responses found in the $i$th partition and $\mathbf{y}_{i}^{\star}=[y^{\star}_{ij}]_{j=1}^{r}$ where $y_{ij}^{\star}$ is the number of observations of $\mathbf{y}_{i}$ that fall into region $j$ (i.e. closest to $z_{ij}$), then the log-likelihood contributions of the $n_{i}$ observations is

[TABLE]

where $\mathbf{f}_{ij}$ is the $j$th element of $\mathbf{f}_{i}$. Consequently, the posterior distribution of $\mathbf{f}_{i}$ can be expressed

[TABLE]

and the marginal distribution after integrating out $\mathbf{f}_{i}$ is

[TABLE]

Given a tessellation, $T$, and using the independence assumption among the partition regions, we can express the marginal likelihood as $p(\mathbf{y}\mid\bm{\theta},[Z_{i}]_{i=1}^{M},T)=\prod_{i=1}^{M}{p(\mathbf{y}_{i}\mid\bm{\theta}_{i},Z_{i})}$. The conditional distribution of the tessellation parameters, $T$, can be obtained by marginalizing $\mathbf{f}_{i}$ in each region:

[TABLE]

This marginalization makes a collapsed MCMC algorithm possible. However, due to lack of conjugacy, $p(\mathbf{y}_{i}\mid\bm{\theta}_{i},Z_{i},T)$ is not available explicitly. Therefore, we use a Laplace approximation (Riihimäki et al., 2014) to integrate over $\mathbf{f}_{i}$ in each partition element. The posterior approximation of $\mathbf{f}_{i}$ takes the following form in each partition

[TABLE]

where $\hat{\mathbf{f}}_{i}=\arg\max_{\mathbf{f}}p(\mathbf{f}_{i}\mid\mathbf{y}_{i},Z_{i},\bm{\theta}_{i})$ and $\Sigma_{i}=((K_{i}+H_{i}BH_{i}^{T})^{-1}+n_{i}(\text{diag}(\mathbf{u}_{i})-\mathbf{u}_{i}\mathbf{u}_{i}^{T}))^{-1}$ with the entries of the vector $\mathbf{u}_{i}$ as $u_{ij}=\exp(\hat{\mathbf{f}}_{ij})/\sum_{j=1}^{r}\exp(\hat{\mathbf{f}}_{ij})$. The Laplace approximation yields an approximate form of the marginal of $p(\mathbf{y}_{i}\mid Z_{i},\bm{\theta}_{i})$ which can be used to approximate $P(T\mid\mathbf{y},\bm{\theta},[Z_{i}]_{i=1}^{M})$ in the MCMC algorithm. A key strategy in Kim et al. (2005) is to marginalize out all parameters which are specific to individual partitions and develop a collapsed MCMC algorithm to obtain posterior samples of the tessellation, $T$. Using the posterior samples of $T$, the best partition (in terms of either posterior probability or marginal likelihood of $\mathbf{y}$) is selected and used to make inference and predictions.

Once a posterior partition is selected, the Laplace approximation makes density estimation within each partition simple and straightforward. Obtaining $\hat{\mathbf{f}}_{i}$ can generally be done using Newton’s method, and is much faster computationally than using MCMC methods to draw from the posterior of $\mathbf{f}_{i}$. Once $\hat{\mathbf{f}}_{i}$ has been found, the posterior density of $p(y_{i})$ is found by simulating from $N(\mathbf{f}_{i}\mid\hat{\mathbf{f}}_{i},\Sigma_{i})$ and applying the log-density transform.

We now discuss the form of $\kappa_{\bm{\theta}_{i}}(\cdot,\cdot)$ and the selection of $\bm{\theta}_{i}$ in each partition. For the purposes of this paper, we assume the covariance function $\kappa_{\bm{\theta}_{i}}(\cdot,\cdot)$, which depends on hyperparameters $\bm{\theta}_{i}=(\sigma_{i}^{2},l_{i})$, is the stationary squared exponential covariance function

[TABLE]

where $\sigma_{i}^{2}$ is the magnitude hyperparameter and $l_{i}$ is a length-scale hyperparameter which together govern the smoothness properties of $f_{i}$. We place a weakly informative half Student-$t$ distribution with one degree of freedom and a variance equal to 10 for the magnitude parameter and the same prior with a variance of 1 for the length-scale hyperparameter (Riihimäki et al., 2014). We then obtain the maximum a posteriori (MAP) estimate of the posterior mode of $\bm{\theta}_{i}$ using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm and set $\bm{\theta}_{i}$ equal to this value. This approach to selecting and fixing $\bm{\theta}_{i}$ can be viewed as an empirical Bayes prior.

2.3 Markov Chain Monte Carlo algorithm

The posterior distribution of the tessellation structure, $T$, does not have an explicit form, thus obtaining posterior draws of the tessellation proceeds using a reversible jump MCMC algorithm (Green, 1995). A reversible jump MCMC algorithm is necessary due to the varying dimension of the parameter space of the tessellation structure. The algorithm to explore this varying dimensional parameter space proceeds as follows:

Initialize the tessellation structure by choosing one observation’s covariate vector as the initial center of the tessellation and choosing a starting weight $\mathbf{w}$ (a simple default is to set $\mathbf{w}=(p^{-1},\ldots,p^{-1})$). 2. 2.

When $1<M<M_{\text{max}}$, with probability 1/4, either add (birth step), delete (death step), or move (moving step) a tessellation center; otherwise, modify the weight vector $\mathbf{w}$. A new center is chosen uniformly from the locations (observed covariate vectors) that are not currently a tessellation center, a center is deleted by randomly choosing an existing center and removing it. Moving a tessellation center entails uniformly choosing one existing center and randomly moving it to another observation’s covariate vector which is not currently a tessellation center. To change the weight vector, a new weight vector $\mathbf{w}^{(p)}$ is proposed by drawing from $q(\mathbf{w}^{(p)}\mid\mathbf{w})\equiv\text{ Dirichlet}(dw_{1},\ldots,dw_{p})$, where $\mathbf{w}=(w_{1},\ldots,w_{p})$ is the vector of the current tessellation weights in the MCMC chain and $d>0$ is a tuning parameter which determines the variance of the proposal distribution $q(\cdot\mid\cdot)$. Accept the proposed change to the tessellation with probability $\alpha$, defined below. 3. 3.

Repeat step 2 until the MCMC chain has converged and enough samples have been generated. Discard burn-in period.

Let $T_{p}$, $T_{c}$ be the proposed and current tessellation structures with $M_{p}$ and $M$ partitions, respectively. We accept the birth, move, death, or change in weight step using the acceptance probability

[TABLE]

It is important to note that for most non-boundary cases, the prior on the tessellation structure $P(T_{p})$ and $P(T_{c})$ does not appear in $\alpha$ due to cancellations with itself and/or the proposal distribution for the birth, death, move, and changing weight steps in the reversible jump MCMC algorithm. When $M$ is at or near the boundary, however, adjustments to $\alpha$ need to be made in order to maintain the reversibility of the MCMC chain. When $M=1$ and a birth step is proposed, we must multiply the ratio in (2) by 3/4. When $M=2$ and we propose a death step, the ratio must be multiplied by 4/3. The reverse must be applied when $M=M_{\text{max}},\text{ }M_{\text{max}}-1$. Note that if the weight vector $\mathbf{w}$ remains unchanged (i.e. birth, death, and move steps), the ratio involving $q(\cdot\mid\cdot)$ will simply be 1.

For a given tessellation drawn from the posterior distribution of $T$, computing the density in each partition of the tessellation is done by simply computing a univariate density estimation using the logistic Gaussian process model. One can then choose the tessellation (partition) with highest posterior probability or marginal likelihood of $y$ and use this as the final density regression model.

3 Some Results on Convergence

In this section, we show the consistency of the proposed conditional density estimation method. If the proposed model is true, then as $n$, the number of points goes to infinity, the posterior density concentrates near a small total variation neighborhood around the true density. We first state the result for the Euclidean norm, denoted by the metric $d$. For a general weighted norm, the result is extended in the appendix.

To prove this result, we first show that the partitions formed by a Voronoi tessellation can adequately approximate the true partition. Then, we show that we have sufficient prior probability for the approximating partition and around any small neighborhood of the true Gaussian process path in supremum norm. Finally, if we have sufficient prior mass around the true density, the likelihood pulls the posterior density towards the data generating density under the true model.

Let ${\mathbf{c}}_{1},{\mathbf{c}}_{2},\dots,{\mathbf{c}}_{m}$ be the centers of some Voronoi tessellation and $R_{1},R_{2},\dots,$ $R_{m}$ be the corresponding Voronoi regions in $\Omega$, a subset of $\mathbb{R}^{d}$ with associated Lebesgue measure $\mathscr{L}$. Let $V_{1},\dots,V_{k}$ be any given partition of $\Omega$. We assume that each region, $V_{i}$, is a finite union of rectangular regions. Our result holds for a general region approximated by a finite union of rectangles. We first prove that any aforementioned region $V_{i}$, can be approximated by the regions of a Voronoi tessellation.

Proposition 1

Given $\epsilon_{1}>0$ there exists $M$ and ${\mathbf{c}}_{1},\dots,{\mathbf{c}}_{M}$ and a partition $J_{1},\dots,J_{k}$ of $\{1,\dots,M\}$ such that $U_{l}=\cup_{i\in J_{l}}R_{i}$ and $\sum_{l=1}^{k}\mathscr{L}(U_{l}\Delta V_{l})\leq\epsilon_{1}$, where $\Delta$ denotes the symmetric differences of sets.

In our proposed method, we use the observed values of the covariates for the centers of the tessellation. Next, we show that a small perturbation of ${\mathbf{c}}_{1},\dots,{\mathbf{c}}_{M}$ from Proposition 1 does not change the partition dramatically and provides an approximation for regions $V_{1},\dots,V_{k}$. Then, we show that any small neighborhood of ${\mathbf{c}}_{1},{\mathbf{c}}_{2},\dots,{\mathbf{c}}_{M}$ contains observed covariates with probability 1 as $n$ goes to infinity. For that we assume that the probability measure on ${\bf x}$, $\tilde{H}(\cdot)$ has a strictly positive, bounded density function. This conclusion implies that we can use covariate points as the centers for the proposed Voronoi tessellation to approximate true partition of the covariate space. We summarize these two results in the form of two following propositions.

Proposition 2

Given $\epsilon_{1}>0$, ${\mathbf{c}}_{1},\dots,{\mathbf{c}}_{M}$ and $R_{1},\dots,R_{M}$ from Proposition 1, we can have $\delta>0$ and Voronoi centers $\mathbf{c}^{\prime}_{1},\dots,\mathbf{c}^{\prime}_{M}$ and corresponding $R^{\prime}_{1},\dots,R^{\prime}_{M}$ such that if $d(\mathbf{c}_{i},\mathbf{c}^{\prime}_{i})<\delta$ then $\sum_{l=1}^{k}\mathscr{L}(U^{\prime}_{l}\Delta V_{l})\leq 2\epsilon_{1}$, where $U^{\prime}_{l}=\cup_{i\in J_{l}}R^{\prime}_{i}$ and $d$ denotes the distance under Euclidean norm.

Proposition 3

Under the setup of Proposition 1, as $n\rightarrow\infty$, we observe $\ \mathbf{x}_{j}\in(\mathbf{c}_{j}\pm\frac{\delta}{2})$ for all $1\leq j\leq M$ with probability 1.

Let

[TABLE]

$\mu_{j}$ is the mean function given in (1) and $f^{*}$ denote the true density. Let $|\eta_{j}^{*}(y)|<k_{0},k_{0}>1$, corresponding to $f^{*}$. We also assume $\eta_{j}^{*}$ to be smooth. Let $\epsilon_{2}<2k_{0}\epsilon_{1}$. We drop the subscript $\mathbf{x}$ from $f$ for notational convenience. Let $V_{1}^{*},\dots,V_{k}^{*}$ be the true underlying partition of $\Omega$ and from Proposition 1 there exists $U_{1}^{*},\dots,U_{k}^{*}$ from the Voronoi approximation. Let $|\mu_{j}^{*}|<k_{0}$ be the true mean function in the $j$th region. Consider the following neighborhood in supremum norm ($\|\|_{\infty}$)

[TABLE]

[TABLE]

Proposition 4

For $\{\mu_{j},\eta_{j}\}$ pairs such that $\mu_{j}\in N_{2}$ and $\eta_{j}\in N_{1}$, $\forall j$; we have $\int|f(y)-f^{*}(y)|<k\epsilon_{2}$, for some $k>0$.

The covariance kernel between points $s$ and $t$ can be written as $K(s,t)=\sigma^{2}K_{0}(\frac{1}{l_{i}}s,\frac{1}{l_{i}}t)$ for the $i$th region $R_{i}$, where $K_{0}(\cdot,\cdot)$ is a smooth kernel. We assume the following,

• $\mathbf{A1}:\log(\text{max }\{\Pi(\sigma>\lambda_{n}),\Pi(\frac{1}{l_{i}}>\nu_{n})\})=O(-n)$.
• $\mathbf{A2}:M_{n}^{2}\lambda_{n}^{-2}\nu_{n}^{-2\alpha}/n\rightarrow\infty$.
• $\mathbf{A3}:{M_{n}}^{1/\alpha}=O(n^{\gamma}),0<\gamma<1$.

Here, $M_{n}$ is of polynomial order of $n$, $\lambda_{n}$ and $\nu_{n}$ are two sequences of constants, and $\alpha\geq 1$, an integer.

Let $\eta()$ be any Gaussian process path, under the smoothness of the covariance kernel the paths are smooth and the derivative process is again a Gaussian process. For any density based on $m\leq M$ partitions we have an $m$-dimensional product function space. We construct sieves on the function space where the probability outside the sieves decreases exponentially with $n$ and establish an entropy bound for the sieves. We use this construction to prove our following convergence result.

Theorem 1

Let $U_{\epsilon^{\prime}}=\{f(y):\int|f(y)-f^{*}(y)|dy<\epsilon^{\prime}\}$. Then, under $\mathbf{A1}-\mathbf{A3}$ and log Gaussian process prior and model (3), $\Pi(U_{\epsilon^{\prime}}\mid\cdot)\rightarrow 1$ with probability one, as $n$ the number of observations goes to infinity.

Even though the main results focus on the neighborhood of the estimated density, the prior favors smaller partitions. Heuristically, if the true partition is further partitioned into smaller partitions, then the true likelihood remains the same over the smaller partitions, but the prior puts $O(n^{-m})$ weight on a partition with $m$ centers. Hence, extra sub-partitions will reduce the posterior probability. Therefore, we should have higher posterior probability for the smaller number of Voronoi centers, as long as it can capture the true data generating partition. We can use a prior satisfying $\mathbf{A1}-\mathbf{A3}$, or truncate the hyperparameters at $\lambda_{n}$ and $\nu_{n}$.

It is important to note that the application of this model differs from the theory in several small ways. First, the applied methodology uses two approximations: the discretized version of $f_{i}$ in each partition, and the Laplace approximation of the marginal of $\mathbf{y}$. The theory, of course, is not based on these practical approximations, but the results must necessarily depend on a reasonable approximation. A measure of the closeness of these approximations to the true underlying model is not undertaken in the present paper, but empirical results of Riihimäki et al. (2014) indicate that these approximations are reasonable in practice for density estimation.

The theory presented provides consistency statements on models which are more general than the applied model. The applied model assumes that the true partition structure is a Voronoi tessellation whereas the theory allows for any partition whose elements are made up of unions of rectangular regions. Modeling this more general structure is beyond the current endeavor of this paper, but opens the door for an even more flexible modeling approach where assuming the true partition structure is a single Voronoi tessellation is too restrictive.

4 Simulations & Applications

4.1 Preliminaries

We implement the partition model on several simulated and real datasets. Since the goal of this methodology is to provide insight as to how the density of $y$ changes with covariates $\mathbf{x}$, the maximum number of partitions for these simulations and applications is capped at 10 to aid in interpretation. Each MCMC chain was run for 10,000 iterations with a burn in period between 1,000 and 2,000 iterations. The data are also centered and scaled to provide numerical stability and allow greater interpretability of the weight vector, $\mathbf{w}$, which performs variable selection on the standardized covariates.

The code to run this algorithm was written in Matlab and utilizes portions of the excellent code developed by Vanhatalo et al. (2013). A basic function to carry out these types of models is available at https://github.com/gitrichhub/bayes-cde.

4.2 One partition

We begin with a very simple simulation with $n=1,000$ data points generated from the following model:

[TABLE]

Note that in this case the simulated data $\mathbf{y}$ have no relationship with $\mathbf{x}$, and therefore the method should favor no splitting of the predictor space (i.e. $M=1$). The MCMC algorithm assigned a 99.85% posterior probability of no splitting of the covariate space and .15% probability of having a partition with two regions, indicating that the model appropriately identified the appropriate relationship between $y$ and $x_{1},\ x_{2}$.

4.3 Piecewise regression

Data ($n=10,000$) was generated by simulating $x_{1}\sim U(0,5),\ x_{2}\sim N(3,2^{2})$, and $y\sim N(f(x_{1}),.25^{2})$ where

[TABLE]

Note in this case we have introduced a covariate, $x_{2}$, which does not have a relationship with $y$. As expected, the weight for variable $x_{2}$ was extremely small (less than $1.5*10^{-11}$) for the entire MCMC chain after burn in, indicating the weight vector works well as a variable selector in this framework.

Figure 2 plots $y$ against $x_{1}$. Colors represent distinct regions from the posterior tessellation with the highest marginal likelihood of $y$. The method successfully identified that the distribution of $y$ is the same for $x<2.5$ and successfully captured changes in the mean function $f(x)$ for $x>2.5$. If the maximum number of allowed partitions is increased to 100, then 65-70 partitions are chosen by the model resulting in a much finer partition over $x>2.5$.

4.4 Bivariate surface

In this simulation we move to the case where the mean of $y$ is a function of two covariates. Data ($n=10,000$) was generated from the following model:

[TABLE]

For a given value of $x_{1}$, the mean function of $y$ is a straight line in $x_{2}$ with slope $[1+e^{-x_{1}}]^{-1}$. Small values of $x_{1}$ correspond to a slope in the $x_{2}$ direction near 0, and large values limit to a slope of 1 in the $x_{2}$ direction.

The contours of the mean function of $y$ are plotted in panel A of Figure 3 and the posterior partition (with highest marginal likelihood of $y$) for $x_{1}$ and $x_{2}$ is plotted in panel B, with colors indicating distinct regions of the tessellation. The tessellation successfully captures the general features of the mean function of $y$. When $x_{1}$ is small, there is little change in the mean of $y$, and the tessellation structure assigns large areas of $\mathbf{x}$ to the same region. As $x_{1}$ increases, we see that a finer partition is induced as the mean structure changes more rapidly over $x_{2}$.

4.5 Changing parametric form

The previous simulated examples have dealt with the case where the distribution of $y$ is related to $x$ only through the mean of $y$. This example allows the mean function, parametric form, and parameters of the distribution of $y$ to vary with $\mathbf{x}$. Specifically,

[TABLE]

where $Z_{x_{2}}\sim\text{Gamma}(2,x_{2})$ with a location shift such that $EZ_{x_{2}}=0$. Thus, for a given value of $x_{1}\geq 5$, the distribution of $y$ is a mixture of Gamma random variables. In this simulation, $x_{1}\sim U(0,10),\ x_{2}\sim U(0,5),\text{ and }x_{3}\sim N(0,5^{2})$. As before, $n=10,000$ and the maximum number of partitions was set to 10. A plot of $y$ against $x_{1}$ is given in Figure 4.

Again, we have a variable, $x_{3}$, which is not related $y$. The weights of the three variables in the posterior partition with highest marginal likelihood of $y$ reflect the fact that $x_{1}$ is the dominant predictor with $x_{2}$ playing a much smaller role and $x_{3}$ playing essentially no role ($\mathbf{w}=(.9814,.0185,.0001)$).

Figure 5 shows the posterior partition with highest marginal probability of $y$. The method successfully identifies the partition for $x_{1}<5$. Interestingly, when the mean of $y$ is changing slowly (when $x_{1}$ is just above 5), the model identifies $x_{2}$ as having an important effect (i.e. for a given $x_{1}$ value between 5 and 7, the partition changes as a function of $x_{2}$). When $x_{1}$ is larger (and the mean of $y$ is changing faster), $x_{2}$’s effect is considered negligible compared to that of $x_{1}$ (i.e. the partition doesn’t change over $x_{2}$ when $x_{1}$ is large).

4.6 Melbourne temperature data

Hyndman et al. (1996) introduced several interesting datasets which are well suited for conditional density estimation methods. One of these is the Melbourne temperature dataset which records each day’s maximum temperature (in Celsius) between 1981 and 1990 in Melbourne, Australia. Interestingly, when each day’s maximum temperature is plotted against the previous day’s maximum temperature, the temperatures fork as the temperature on the x-axis increases, making it unsuitable for usual regression or other non-parametric regression techniques which generally assume a normally distributed error.

The partition model was applied to this dataset, and the most common partition is shown in Figure 6. Figure 7 shows the estimated density of $y$ in several of the partitions. These densities reveal a steady progression towards bimodality in the distribution of today’s temperature as yesterday’s temperature increases. Hyndman et al. (1996) explains that this bimodality is a consequence of high pressure systems passing over the city, which are sometimes followed by cold fronts during hotter months (resulting in significantly cooler days following hot days). Thus hot days are followed by hot days if there is no cold front or significantly cooler days when a cold front passes over Melbourne.

4.7 Dow Jones changepoint analysis

The partition model framework can also be used to perform changepoint analyses. James and Matteson (2013) perform a changepoint analysis on the weekly log returns of the Dow Jones Industrial Average (DJIA) index. Applying our partition model framework to this dataset yields 6 partitions (5 changepoints). The returns are plotted in Figure 8 with colors representing the posterior tessellation with highest marginal likelihood. James and Matteson (2013) identify 4 changepoints in their analysis, of which 3 correspond closely with the results from our model (12/9/1996, 4/14/2003, and 10/22/2007). The other two changepoints from our model are 2/25/1991 and 5/27/2002.

The selected changepoints have some interesting relationships with world market conditions. The changepoint in February 1991 corresponds to the Japanese asset price bubble. The changepoint identified in May of 2002 corresponds to the U.S. stock market downturn of 2002. The changepoint in April of 2003 occurs during the U.S. invasion of Iraq, and the October 2007 changepoint is near the time of sub-prime mortgage crisis. The right panel of Figure 8 shows the posterior mean densities of the returns in each partition. There appears to be three groups with similar periods of volatility: February 1991 to December 1996 and April 2003 to October 2007 appear to be marked by lower levels of volatility whereas the period between May 2002 and April 2003 is marked by higher volatility. The remaining time periods appear to have similar levels of volatility.

One major advantage of this method in changepoint analysis is its ability to estimate the density in each region of the partition. We are not limited to assuming any parametric form of the density in each region, and therefore our the estimates of the density of the returns are not restricted to Gaussian or symmetric heavy tailed distributions. Indeed, the period between 5/27/2002 and 4/14/2003 is slightly left skewed.

4.8 Windmill data

Conditional density estimation is particularly useful when it is unclear how the density of $y$ changes with respect to the predictors, $\mathbf{x}$. The exact relationship of electrical power output in windmills with various covariates (wind speed, wind direction, air density, wind sheer, & turbulence intensity) is unknown, and does not appear to follow any known parametric form.

We analyze a wind turbine dataset in order to predict the average power output given a set of predictors. Using a random subsample of 10,000 observations from a larger wind turbine dataset, the partition model was fit to the data using 5 covariates (wind speed, wind direction, air density, wind sheer, & turbulence intensity). Density estimates were fit in each region of the tessellation with the highest marginal likelihood.

Figure 9 shows a plot of the posterior partition structure with two of the covariates, wind speed and direction. From the plot, it is easy to see the importance of wind speed and direction in determining the overall partition structure. The overlap in regions in this two-dimensional view of the partition indicates the effect of other covariates in determining the 5-dimensional partition structure.

Figure 10 shows posterior densities (with 90% confidence bands) in four regions of the tessellation. Note how the distribution of power output changes dramatically throughout the covariate space, demonstrating the need for a density regression technique to more accurately determine the density of $y$ in various regions of $\mathbf{x}$.

5 Conclusion

The combination of logistic Gaussian process density estimation with partition modeling provides an excellent framework for determining how and where the density of $y$ changes over the covariate space. Furthermore, it has desirable consistency properties. The applications of this model will help analysts determine which variables are important in predicting the density of $y$, as well as better estimate the density of $y$ over the covariate space. This flexible method will help give greater insight into datasets in which the relationship between the density of $y$ and $\mathbf{x}$ is unknown or difficult to formulate parametrically. As part of our contribution, we provide the first (to our knowledge) publicly available software to perform Bayesian conditional density estimation. The code is publicly available at https://github.com/gitrichhub/bayes-cde.

6 Appendix

We show the result from Proposition 1 and 2 for any rectangular region in $d$-dimensional covariate space. We show that any rectangular region can be approximated up to any accuracy by a union of Voronoi regions. The result is shown using a grid-based tessellation construction. For general regions, we extend the result thereafter. We then prove our main theorem and address the Voronoi region construction issue under the weighted norm. Later, we address a special case where the true partition in the covariate space is a Voronoi tessellation.

6.1 Proof of Proposition 1

Let $\xi>0$ be any generic small constant. Let $R$ be any $d$-dimensional rectangle with minimum edge length $l$ and maximum edge length $L$. Assume $l>8d\xi$. We construct two outer rectangles $R^{(1)},R^{(2)}$ with edge length $L_{i}+4d\xi$ and $L_{i}+8d\xi$, and two inner rectangles $R_{(1)},R_{(2)}$ with edge lengths $L_{i}-4d\xi$ and $L_{i}-8d\xi$, where $L_{i}$ is the length of $i$th edge.

The following gives the explicit description of $R^{(1)}$ and $R_{(1)}$. Rectangles $R^{(2)}$ and $R_{(2)}$ can be constructed similarly. In particular, let $\overline{A_{i}B_{i}}$ be the edge with length $L_{i}$ and $O$ be the center of mass of the rectangle. Let, $\overline{A_{i}^{\prime}B_{i}^{\prime}}$ be the edge of the outer rectangle with same center of mass $O$ and length $L_{i}+4d\xi$, and $\overline{a_{i}^{\prime}b_{i}^{\prime}}$ is the edge of inner rectangle with center of mass $O$ and length $L_{i}-4d\xi$. Here, $A_{i}^{\prime}$, $a_{i}^{\prime}$ lie on the half-line $\overrightarrow{OA_{i}}$, and $B_{i}^{\prime}$, $b_{i}^{\prime}$ lie on the half-line $\overrightarrow{OB_{i}}$, and $\overleftrightarrow{A_{i}^{\prime}B_{i}^{\prime}}$, $\overleftrightarrow{A_{i}B_{i}}$ and $\overleftrightarrow{a_{i}^{\prime}b_{i}^{\prime}}$ are parallel to each other.

For choosing Voronoi centers, we construct $\xi$-distanced equi-spaced grids at each edge of $R$ and use the induced grid points in the interior of $R$ associated with $d$-dimensional rectangles with edge length $\xi$. The number of centers that fall in the rectangle $R$ with edge length $L_{i}$ is of the order of $(L/\xi)^{d}$. Let, ${\bf c}_{1},\dots,{\bf c}_{M_{0}}$ be the centers that fall within $R_{(1)}$. Let $R_{1},\dots,R_{M_{0}}$, be the corresponding Voronoi regions.

The proposed Voronoi tessellation has centers at ${\bf c}_{1},\dots,{\bf c}_{M_{0}}$ and at grid-centers that lie on the sides of $R$. Note that for any point in $R_{(1)}$, the corresponding Voronoi center cannot be outside $R_{(1)}$. For any such point, the Voronoi centers that lie on the sides of $R$ have at least $2\xi d$ distance. At the same time, for any point in $R_{(1)}$, that point lies in a $\xi$ edge length $d$-dimensional rectangular grid, where its vertices are in the set $\{{\bf c}_{1},\dots,{\bf c}_{M_{0}}\}$. Hence, we have a vertex at maximum $\xi d$ distance (a very conservative bound). Let, $\{{\bf C}\}=\{{\bf c}_{1},\dots,{\bf c}_{M_{0}}\}$ and let $\{\tilde{\bf C}\}$ be the union of $\{{\bf C}\}$ and the grid-centers that lie on the sides of $R$.

Similarly, for any point outside $R$, the corresponding Voronoi center has to be outside of $R_{(1)}$. We use regions corresponding to $\{{\bf C}\}$ to approximate $R$.

Note $\cup_{i=1}^{M_{0}}R_{i}=R^{*}\subset R$ and $\mathscr{L}(R^{*}\Delta R)\leq\mathscr{L}(R\Delta R_{(1)})\leq 8d\xi L$. Here, $\xi$ can be arbitrarily small, thus proving our claim.

6.2 Proof of Proposition 2

From the proof of Proposition 1, we shift each center $\bf{c}_{1},\dots,\bf{c}_{M_{0}}$ by less than or equal to $\delta$ distance in Euclidean norm.

After the perturbation let ${\bf c}_{i}^{*}$, $i=1,\dots,M_{0}$ be the new Voronoi centers. Then, for $\delta<.5\xi$, and for any point in $R_{(1)}$, the corresponding Voronoi centers must be in ${\bf c}^{*}_{1},\dots,{\bf c}^{*}_{M_{0}}$ as the distance between cannot be more than $d\xi+.5\xi$, which follows from the following argument.

As earlier in Proposition 1, points in $R_{(1)}$ will be inside an $\xi$ edge length d-dimensional rectangle induced by the grid construction. Therefore, points in $R_{(1)}$ will have a vertex from the set $\{{\bf c}_{1},\dots,{\bf c}_{M_{0}}\}$ within $\xi d$ distance. Also, there exists ${\bf c}^{*}_{i}$’s within $.5\xi$ distance from the grid vertices in $R_{(1)}$. Hence, from the triangle inequality, for any point in $R_{(1)}$, there exists a point in $\{{\bf c}^{*}_{1},\dots,{\bf c}^{*}_{M_{0}}\}$ within distance $d\xi+.5\xi$.

From any grid point on the sides of $R$, the distance of any point in $R_{(1)}$ is at least $2d\xi$ and hence, after perturbing these grid points by $\delta$, we have a distance of at least $2d\xi-\delta>d\xi+.5\xi$ as $d\geq 1$ and $\delta<.5\xi$. Therefore, points in $R_{(1)}$ will be in Voronoi regions corresponding to the centers ${\bf c}^{*}_{1},\dots,{\bf c}^{*}_{M_{0}}$.

Similarly, for any points outside $R^{(1)}$, the Voronoi centers cannot be one of $\bf{c}^{*}_{1},\dots,\bf{c}^{*}_{M_{0}}$. As, $\mathscr{L}(R^{(1)}\Delta R_{(1)})$ $\leq 16d\xi L$, it concludes our claim.

6.3 General regions

We have shown the result for a rectangular region, that is, the covariate space is partitioned into a rectangle and its outside. For general regions, the proof follows from writing each partition as a union of non-intersecting finitely many rectangles.

Suppose, each partition of the covariate space can be written as a non-intersecting union of finitely many rectangles. Let $R^{*}_{1},\dots,R^{*}_{M}$ be the $M$ generating rectangles. Let $\{\tilde{\bf C}_{i}\}$ be the Voronoi centers constructed for $R^{*}_{i}$, as in the proof of Proposition 1 (using the grid points on the interior rectangle and the grid points that lie on the side of the given rectangle $R_{i}^{*}$). Then, $\cup_{i}\{\tilde{\bf C}_{i}\}$ gives the Voronoi tessellation for the covariate space, where each partition element can be approximated by a union of Voronoi regions, with approximation error $O(\xi)$ as in Proposition 1.

For unions of rectangles, we use the union of regions corresponding to Voronoi centers $\cup_{i\in I}\{{\bf C}_{i}\}$ to approximate $\cup_{i\in{\bf I}}R^{*}_{i}$, for any index set ${\bf I}$, where the approximation error for each $R^{*}_{i}$ is $O(\xi)$. Hence, the approximation error is of the order $O(\xi)$ for $\cup_{i\in{\bf I}}R^{*}_{i}$ , for any ${\bf I}$. Therefore, the Lebesgue measure of the Voronoi region with respect to ${\bf B}_{c}=\cup_{i\in I}\{\tilde{\bf C}_{i}\}\backslash\cup_{i\in I}\{{\bf C}_{i}\}$, that is the centers on the boundaries of the rectangles, is $O(\xi)$. Hence, if $\Omega=\cup_{i=1}^{k}V_{i}$, where, $V_{l}=\cup_{i\in{\bf I}_{l}}R^{*}_{i}$ for disjoint partitions of ${I}_{l}$ of $\{1,\dots,k\}$. Then regions with respect to $\cup_{i\in I_{l}}\{{\bf C}_{i}\}$’s are used to approximate $V_{l}$’s; $l<k$ and regions corresponding to $\cup_{i\in I_{k}}\{{\bf C}_{i}\}\cup{\bf B}_{c}$ approximate $V_{k}$.

6.4 Proof of Proposition 3

This claim follows as the covariate distribution has positive density over the underlying domain. Therefore, relative frequency of the observations in any open set converges almost surely to the corresponding measure of the open set under measure $\tilde{H}(\cdot)$. Therefore, each small neighborhood around ${\bf c_{j}}$ has some observation from the covariate vector with probability one, as $n$ the number of observations, goes to infinity.

In particular let $\tilde{O}$ be any open set and suppose $\tilde{H}(\tilde{O})=\xi^{\prime}>0$, where $\tilde{H}$ is the underlying measure for the $d$-dimensional covariate distribution. Then, $p_{i}=P(\text{ there is no observation in }\tilde{O}\text{ for }{\bf x}_{i};\ i=1,\dots,n)=(1-\xi^{\prime})^{n}$.

As, $\sum p_{i}<\infty$, then by the Borel-Cantelli lemma, with probability one, there are observations in $\tilde{O}$ for all but finitely many $n$’s.

6.5 Proof of Theorem 1

6.5.1 Sieve Construction

Let $\eta(\cdot)$ be any Gaussian process path and $\eta_{j}(\cdot)$ be corresponding to the $j$th Voronoi region. For any $j$th Voronoi region, let $\mathbf{K}^{n}_{j}$=$\{\|D^{l}(\eta_{j}(.))\|_{\infty}<M_{n},l\leq\alpha,\sigma_{i}<\lambda_{n},\frac{1}{l_{i}}<\nu_{n}\}$ and $\mathbf{K}_{\beta_{j}}^{n}=\mathbf{K}_{b}^{n}=\{\|\mathbf{\beta}_{j}\|_{\infty}<\sqrt{n}\}$, where $D^{l}$ implies the $l$th derivative. Then the corresponding log covering number of smooth Gaussian process paths under supremum norm is $log(N(\epsilon,\mathbf{K}_{j}^{n},\|\|_{\infty}))=o(n),$ and $log(\Pi((\mathbf{K}^{n}_{j})^{c}))\leq O(-n)$ and $log(\Pi((\mathbf{K}_{b})^{c}))=O(-n)$ (see Tokdar and Ghosh (2007); Ghosal and Roy (2006), Lemma 2).

Also, covering $\mathbf{K}_{b}^{n}$ with an $\epsilon$ space grid, the covering number is a polynomial of $n$ of order 1, and $log(N(\epsilon,\mathbf{K}_{b}^{n},\|\|_{\infty}))=o(n)$.

Let, $\mathbf{K}^{n}_{j,\beta}$ be the prior product space of $\mathbf{K}_{\beta_{j}}^{n}$ and $\mathbf{K}^{n}_{j}$. Covering the product space by the product rectangular set $\epsilon$-width grids of $\mathbf{K}_{\beta}^{n}$ and $\mathbf{K}^{n}_{j}$, from the corresponding spaces, we have the covering number of the product space as the product of two covering numbers. Hence, we have

[TABLE]

Also, $log(\Pi((\mathbf{K}^{n}_{j,\beta})^{c}))=O(-n)$.

6.5.2 Entropy bound for the Voronoi tessellation

If for the $j$th tessellation, for the coefficients $\|\mathbf{\beta}_{1,j}-\mathbf{\beta}_{2,j}\|_{\infty}<\epsilon$, then for the corresponding mean functions, $\|\mu_{j,1}-\mu_{j,2}\|_{\infty}<\gamma_{1}\epsilon$ for some $\gamma_{1}>0$, as $y$’s are supported on bounded regions. Hence, if $\|\mathbf{\beta}_{1,j}-\mathbf{\beta}_{2,j}\|_{\infty}<\epsilon$ and $\|\eta_{j,1}-\eta_{j,2}\|_{\infty}<\epsilon$, then for the corresponding densities $f_{1}$ and $f_{2}$, we have $KL(f_{1},f_{2})<\gamma_{2}\epsilon$, $\int|f_{1}-f_{2}|dy<\gamma_{2}\epsilon$; $\gamma_{2}>0$ (Proposition 4). Choosing, $\alpha\epsilon$ instead of $\epsilon$, with sufficiently small $\alpha>0$, we have $\int|f_{1}-f_{2}|dy<\epsilon$. Hence, writing in terms of $L_{1}$ covering number,

[TABLE]

Let, $\mathscr{K}_{m}$, $m<M_{max}$, be the m-dimensional product space for some Voronoi tessellation, $\mathscr{K}_{m}=\prod_{l=1}^{m}\mathbf{K}^{n}_{l,\beta}$, corresponding to a Voronoi tessellation based on $m$ centers. Let $\mathscr{K}$ be the union of all such ${n\choose m},\ m\leq M_{max}$, many combinations of $\mathscr{K}_{m}$’s for all possible Voronoi center selection. As, $log({n\choose m})=o(n)$, then using equation (5), we have:

[TABLE]

where $N(\epsilon,\mathscr{K},\|\|_{1})$ is the $L_{1}$ covering number of $\mathscr{K}$.

6.5.3 Combining the parts

Let $U_{\epsilon^{\prime}}$ be an $\epsilon^{\prime}$ radius $L_{1}$ ball around $f^{*}$. Then,

[TABLE]

Here, $\Pi$ denotes the prior distribution and ${\bf y}$ denotes the $n$ length observation vector. Let, $N_{1,2}^{\epsilon_{2}}$ be the set $\eta_{j}\in N_{1}$ and $\mu_{j}\in N_{2}$, from Proposition 4. For any fixed $\alpha_{1},\epsilon_{3}>0$, choosing $\epsilon_{2}$ small enough, for $\eta_{j}()$ and $\bf{\beta}_{j}$ in $N_{1}$ and $N_{2}$, and for the corresponding density $f$, we have $KL(f^{*},f)<\alpha_{1}\epsilon_{3}$. Also, $\Pi(N_{1,2}^{\epsilon_{2}})>0$. (Tokdar and Ghosh (2007))

This is similar to settings of Tokdar and Ghosh (2007), other than an extra prior term $e^{-Mlogn}$ in the denominator corresponding to the Voronoi center selection probability, which goes to zero at a polynomial rate. Note that, $e^{n\alpha_{1}\epsilon_{3}}e^{-Mlogn}>1$ as $n$ goes to infinity for any $\alpha_{1}>0$. Hence, the proof follows. We give a brief sketch in the following argument. See Ghosh and Ramamoorthi (2003) (Theorem 4.4.3 and 4.4.4 proof) for the details and related test construction that we use next.

Writing,

[TABLE]

where the subscript ${\mathscr{K}}$ denotes the integral over $U_{\epsilon^{\prime}}^{c}\cap{\mathscr{K}}$, and subscript ${\mathscr{K}}^{c}$ denotes the integral over $U_{\epsilon}^{\prime c}\cap{\mathscr{K}^{c}}$.

We can construct test functions $0\leq\Psi_{n}\leq 1$ such that $E_{f^{*}}(\Psi_{n})<e^{-c_{1}n{\epsilon^{\prime}}^{2}/8}$ and $\sup_{f\in U_{\epsilon^{\prime}}^{c}}E_{f}(1-\Psi_{n})<e^{-c_{1}n{\epsilon^{\prime}}^{2}/8}$, for any $c_{1}<1$, for large $n$. We write,

[TABLE]

Then, $E_{f^{*}}((1-\Psi_{n}){{\bf N}^{n}_{1}}_{\mathscr{K}})\leq\sup_{f\in U_{\epsilon^{\prime}}^{c}}E_{f}(1-\Psi_{n})<e^{-c_{1}n{\epsilon^{\prime}}^{2}/8}$ and $E_{f^{*}}[(1-\Psi_{n}){{\bf N}^{n}_{1}}_{\mathscr{K}^{c}}]\leq E_{f^{*}}({{\bf N}^{n}_{1}}_{\mathscr{K}^{c}})<e^{-c_{2}n}$, $c_{2}>0$.

We have, $\Pi(N_{1,2}^{\epsilon_{2}})>0$ and in $N_{1,2}^{\epsilon_{2}}$ the integrand in ${{\bf D}^{n}_{1}}$ is greater than $e^{-2n\alpha_{1}\epsilon_{3}}$ with probability one. Hence, choosing $\epsilon_{3}={\epsilon^{\prime}}^{2}$, $2\alpha_{1}<c_{1}/16$,

[TABLE]

Using Markov inequality and Borel-Cantelli Lemma, as $P(\Pi(U_{\epsilon^{\prime}}^{c}\mid\cdot)>e^{-c_{1}n{\epsilon^{\prime}}^{2}/32}$ \text{infinitely often})$$\leq\lim_{k\rightarrow\infty}\sum_{k}^{\infty}2e^{-c_{1}n{\epsilon^{\prime}}^{2}/32}\rightarrow 0 and $\Pi(U_{\epsilon^{\prime}}^{c}\mid\cdot)$ goes to zero almost surely.

6.6 General weight function for the norm

6.6.1 Voronoi centers derivation

Let $\bf{w}^{*}$ be the true $d$-dimensional weight vector with entries $w_{i}^{*}>0$. Let ${\bf x}_{w}^{*}=T_{w^{*}}({\bf x})$ be the transformed $d$-dimensional scaled covariate vector with $x_{i,w}^{*}={{w}_{i}^{*}}^{\frac{1}{2}}x_{i}$. We can create the Voronoi tessellation given in Proposition 1 and extend it to a general region of a transformed covariate space ${\bf x}_{w}^{*}$ and create the partition, as Euclidean distance in ${\bf x}^{*}$ is equivalent to the weighted metric. Therefore, centers in the covariate space can be achieved by performing an inverse transformation on the centers achieved by this scaling.

Let $\delta_{w}$ be a small constant such that for small $\delta_{w}$ neighborhood of ${\bf w}^{*}$, that is $\|{\bf w}-{\bf w}^{*}\|_{\infty}<\delta_{w},w_{i}^{*}>\delta_{w}\ \forall i$. Let, ${\bf x}_{1}$ and ${\bf x}_{2}$ be two points in the covariate space, and ${\bf x}^{*}_{1}$ and ${\bf x}^{*}_{2}$, and ${\bf x}_{1,w}$ and ${\bf x}_{2,w}$ be their transformed versions for weights ${\bf w}^{*}$ and ${\bf w}$.

Let $d_{w}({\bf x}_{1},{\bf x}_{2})$ be the metric for weight ${\bf w}$. Choosing $\delta_{w}$ small enough, we can have $|d_{w}({\bf x}_{1},{\bf x}_{2})-d_{w^{*}}({\bf x}_{1},{\bf x}_{2})|<\delta$ if $|w_{i}-w_{i}^{*}|<\delta_{w},\forall i$.

Choosing $\delta_{w}$ small, we can have for ${\bf A}\subset\Omega$ : $\mathscr{L}({\bf A})<\epsilon\implies\mathscr{L}(T_{w}({\bf A}))<k_{1}\epsilon$, and $\mathscr{L}(T_{w}({\bf A}))<\epsilon\implies\mathscr{L}({\bf A})<k_{1}\epsilon$; for some constant $k_{1}>0$ for the set $\{{\bf w}:\|{\bf w}-{\bf w}*\|_{\infty}<\delta_{w}\}$. This step follows from a change of variable argument using the fact $w_{i}$’s are uniformly bounded away from zero on a $\delta_{w}$ neighborhood of ${\bf w}^{*}$.

Let ${\bf c}_{1}^{w^{*}},\dots,{\bf c}_{M}^{w^{*}}$ be the Voronoi centers in the covariate space under metric $d_{w^{*}}$ and with regions $R_{1},\dots,R_{M}$ for the setting of Proposition 1. Therefore, under metric $d_{w^{*}}$, for true ${\bf w}^{*}$, we have unions of regions, $U_{l}$’s, approximating the true partition $V_{l}$’s.

Let, $S_{\delta,w}=\{{\bf x}\in\Omega:|d_{w}({\bf x},{\bf c}_{i})-d_{w}({\bf x},{\bf c}_{j})|<2\delta\text{ for some }i,j\}$. Letting $\delta$ decrease to zero, $S_{\delta,w}$ decreases to a set formed by $M\choose 2$ many $d-1$ dimensional hyperplanes which have equal distances from any two Voronoi centers. Each of the hyperplanes have measure zero. Hence, for any $\delta_{1}>0$, we can choose $\delta$ small enough such that $\mathscr{L}(S_{\delta,w^{*}})<\delta_{1}$.

Suppose for any ${\bf w}$ in the $\delta_{w}$ neighborhood of ${\bf w}^{*}$, we have ${\bf x}\in\Omega$ in the Voronoi region corresponding to ${\bf c}_{i}^{w^{*}}$ under $d_{w^{*}}$. Then we show that if ${\bf x}\notin S_{\delta,w^{*}}$, ${\bf x}$ lies in the Voronoi region with the same center under metric $d_{w}$ and vice versa.

The last step follows from the fact that for any point in the covariate space the distances from the centers under two metrics $d_{w}$ and $d_{w^{*}}$ differ by at most $\delta$; hence, if $|d_{w^{*}}({\bf x},{\bf c}_{i})-d_{w^{*}}({\bf x},{\bf c}_{j})|>2\delta$ and $d_{w^{*}}({\bf x},{\bf c}_{i})<d_{w^{*}}({\bf x},{\bf c}_{j})$, then $d_{w}({\bf x},{\bf c}_{i})<d_{w}({\bf x},{\bf c}_{j})$.

Let $\xi>0$ be any small constant. For the centers ${\bf c}_{1}^{w^{*}},\dots,{\bf c}_{M}^{w^{*}}$, we have $R_{1},\dots,R_{M}$ the Voronoi region under $d_{w^{*}}$. Then, choosing $\delta_{w}$ small, for $\tilde{R}_{1},\dots,\tilde{R}_{M}$ the regions under $d_{w}$, we have $\sum R_{i}\Delta\tilde{R}_{i}<\xi$, as we can make the measure of the set $S_{\delta,w^{*}}$ arbitrary small by choosing small $\delta_{w}$.

Suppose, we have $U_{l}$ and $\tilde{U}_{l}$’s corresponding to $d_{w^{*}}$ and $d_{w}$, respectively. Then, $\sum V_{i}\Delta\tilde{U}_{i}<\sum V_{i}\Delta{U}_{i}+\xi$. Hence, choosing $\xi$ small enough, we can use the ${\bf w}$ in the $\delta_{w}$ neighborhood of the true weight ${\bf w^{*}}$ to approximate the covariate space partition using Voronoi regions.

6.6.2 Prior mass condition

Hence, we can use the Voronoi centers for ${\bf w}^{*}$ for a $\delta_{w}$ supremum neighborhood around ${\bf w}^{*}$ for approximating the true partition of the parameter space. Under Dirichlet a prior, that neighborhood has positive probability and hence the proof of Theorem 1 holds.

6.7 Special case: true partition is a Voronoi tessellation

Suppose the true partition of the covariate space is given by $M$ Voronoi regions corresponding to centers ${\bf c}_{1}^{*},\dots,{\bf c}_{M}^{*}$. Let, ${\bf c}_{1},\dots,{\bf c}_{M}$ be in the neighborhood of ${\bf c}_{i}^{*}$ with $d({\bf c}_{i},{\bf c}_{i}^{*})<\delta_{c}$. Then, for any ${\bf x}\in\Omega$, $|d({\bf x},{\bf c}_{i})-d({\bf x},{\bf c}_{i}^{*})|<\delta_{c}$.

Let $S_{\delta_{c}}=\{{\bf x}\in\Omega:|d({\bf x},{\bf c}_{i}^{*})-d({\bf x},{\bf c}_{j}^{*})|<2\delta_{c}\text{ for some }i,j\}$. We choose $\delta_{1}>0$ to be any small constant. Choosing $\delta_{c}$ small enough, we can have $\mathscr{L}(S_{\delta_{c}})<\delta_{1}$ from the argument given in Section 7.6.1 for the general weight function. Outside $S_{\delta_{c}}$, we have $d({\bf x},{\bf c}_{i}^{*})<d({\bf x},{\bf c}_{j}^{*})\implies d({\bf x},{\bf c}_{i})<d({\bf x},{\bf c}_{j})$ and vice versa (triangle inequality). Therefore, for $R_{1},\dots,R_{M}$ the regions corresponding to ${\bf c}_{1}^{*},\dots,{\bf c}_{M}^{*}$ and $\tilde{R}_{1},\dots,\tilde{R}_{M}$ the regions corresponding to ${\bf c}_{1},\dots,{\bf c}_{M}$, we have $\sum R_{i}\Delta\tilde{R}_{i}<\delta_{1}$.

Hence, by picking centers at $\delta_{c}$ neighborhoods around the true Voronoi centers, we can approximate the region with $\delta_{1}$ accuracy. From Proposition 3, we have observations in $\delta_{c}$ neighborhood with probability 1 and can choose our Voronoi tessellation.

6.8 Proof of Proposition 4

This proof follows using techniques similar to Tokdar and Ghosh (2007). A sketch can be given as follows. The set where the two Gaussian process paths are more than $\epsilon$ away in $N_{1}$ has Lebesgue measure of order $\epsilon$ (by the construction of $N_{1}$ with $\epsilon_{2}=\epsilon_{1}=\epsilon$) and hence, its probability under $\tilde{H}$ is less than $k_{3}\epsilon$ where $k_{3}>0$ is fixed (depends on the bounded density function $\tilde{h}$ of $\tilde{H}(.)$). Consider $f_{1}$ and $f_{2}$ with corresponding ${\bf\beta}_{1}$, ${\bf\beta}_{2}$, $\eta_{1}(\cdot)$ and $\eta_{2}(\cdot)$ in $N_{1}$ and $N_{2}$, respectively. Then, for the distribution corresponding to $f_{1}$ and $f_{2}$, $|g_{\bf x}^{f_{1}}(y)-g_{\bf x}^{f_{2}}(y)|=g_{\bf x}^{f_{1}}(y)O(\epsilon)$ outside of an $O(\epsilon)$ measure set of the covariate space. Then, from the fact that $y$ is supported on a bounded region, $|\int g_{\bf x}^{f_{1}}(y)d(x,y)-\int g_{\bf x}^{f_{2}}(y)d(x,y)|=O(\epsilon)$ and the proof follows.