Mass-shifting phenomenon of truncated multivariate normal priors

Let $p(\cdot)$ denote the density of a $\mathcal{N}_{N}(\mathbf{0},\Sigma)$ distribution truncated to the non-negative orthant in $\mathbb{R}^{N}$,

The density $p$ is clearly unimodal with its mode at the origin. However, for certain classes of non-diagonal $\Sigma$, we surprisingly observe that the lower-dimensional marginal distributions increasingly shift mass away from the origin as $N$ increases. This observation is quantified in Theorem 2, where we provide non-asymptotic estimates for marginal probabilities of events of the form $\{\theta_{1}\leq\delta\}$, for $\delta>0$. En-route to the proof, we derive a novel Gaussian comparison inequality in Lemma 1. An immediate implication of this mass-shifting phenomenon is that corner regions of the support $\mathcal{C}$, where a subset of the coordinates take values close to zero, increasingly become low-probability regions under $p(\cdot)$ as dimension increases. From a statistical perspective, this helps explain a paradoxical behavior in Bayesian constrained regression empirically observed in Neelon and Dunson, (2004) and Curtis and Ghosh, (2011), where truncated normal priors led to biased posterior inference when the underlying function had flat regions.

A common approach towards Bayesian constrained regression expands the function in a flexible basis which facilitates representation of the functional constraints in terms of simple constraints on the coefficient space, and then specifies a prior distribution on the coefficients obeying the said constraints. In this context, the multivariate normal distribution subject to linear constraints arises as a natural conjugate prior in Gaussian models and beyond. Various basis, such as Bernstein polynomials (Curtis and Ghosh,, 2011), regression splines (Cai and Dunson,, 2007; Meyer et al.,, 2011), penalized spines (Brezger and Steiner,, 2008), cumulative distribution functions (Bornkamp and Ickstadt,, 2009), restricted splines (Shively et al.,, 2011), and compactly supported basis (Maatouk and Bay,, 2017) have been employed in the literature. For numerical illustrations in this article, we shall use the formulation of Maatouk and Bay, (2017) where various restrictions such as boundedness, monotonicity, convexity, etc were equivalently translated into non-negativity constraints on the coefficients under an appropriate basis expansion. They used a truncated normal prior as in (1.1) on the coefficients, with $\Sigma$ induced from a parent Gaussian process on the regression function; see Appendix A for more details.

To motivate our theoretical investigations, the two panels in Figure 1 depict the estimation of two different monotone smooth functions on $[0,1]$ based on 100 samples using the basis of Maatouk and Bay, (2017) and a joint prior $p(\cdot)$ as in (1.1) on the $N=50$ dimensional basis coefficients. The same prior scale matrix $\Sigma$ was employed across the two settings; the specifics are deferred to Section 3 and Appendix A. Observe that the function in the left panel is strictly monotone, while the one on the right panel is relatively flat over a region. While the point estimate (posterior mean) as well as the credible intervals look reasonable for the function in the left panel, the situation is significantly worse for the function in the right panel. The posterior mean incurs a large bias, and the pointwise 95$\%$ credible intervals fail to capture the true function for a substantial part of the input domain, suggesting that the entire posterior distribution is biased away from the truth. This behavior is perplexing; we are fitting a well-specified model with a prior that has full support¹¹1the prior probability assigned to arbitrarily small Kullback–Leibler neighborhoods of any point is positive. on the parameter space, which under mild conditions implies good first-order asymptotic properties (Ghosal et al.,, 2000) such as posterior consistency. However, the finite sample behavior of the posterior under the second scenario clearly suggests otherwise.

Functions with flat regions as in the right panel of Figure 1 routinely appear in many applications; for example, dose-response curves are assumed to be non-decreasing with the possibility that the dose-response relationship is flat over certain regions (Neelon and Dunson,, 2004). A similar biased behavior of the posterior for such functions under truncated normal priors was observed by Neelon and Dunson, (2004) while using a piecewise linear model, and also by Curtis and Ghosh, (2011) under a Bernstein polynomial basis. However, a clear explanation behind such behavior as well as the extent to which it is prevalent has been missing in the literature, and the mass-shifting phenomenon alluded before offers a clarification. Under the basis of Maatouk and Bay, (2017), a subset of the basis coefficients are required to shrink close to zero to accurately approximate functions with such flat regions. However, the truncated normal posterior pushes mass away from such corner regions, leading to the bias. Importantly, our theory also suggests that the problem would not disappear and would rather get accentuated in the large sample scenario if one follows standard practice of scaling up the number of basis functions with increasing sample size, since the mass-shifting gets more pronounced with increasing dimension. To illustrate this point, Figure 2 shows the estimation of the same function in the right panel of Figure 1, now based on 500 samples and $N=50$ and $N=250$ basis functions in the left and right panel respectively. Increasing the number of basis functions indeed results in a noticeable increase in the bias as clearly seen from the insets which zoom into two disjoint regions of the covariate domain. A similar story holds for the basis of Neelon and Dunson, (2004) and Curtis and Ghosh, (2011).

Neelon and Dunson, (2004) and Curtis and Ghosh, (2011) both used point-mass mixture priors as remedy, which is a natural choice under a non-decreasing constraint. However, their introduction becomes somewhat cumbersome under the non-negativity constraint in (1.1). As a simple remedy, we suggest introducing a multiplicative scale parameter for each coordinate a priori and further equipping it with a prior mixing distribution which has positive density at the origin and heavy tails; a default candidate is the half-Cauchy density (Carvalho et al.,, 2010). The resulting prior shrinks much more aggressively towards the origin, and we empirically illustrate its superior performance over the truncated normal prior. This empirical exercise provides further support to our argument.

In this section, we connect the theoretical findings in the previous section to posterior inference in Bayesian constrained regression models. We work under the setup of a usual Gaussian regression model,

where we assume $x_{i}\in[0,1]$ for simplicity. We are interested in the situation when the regression function $f$ is constrained to lie in some space $\mathcal{C}_{f}$ which is a subset of the space of all continuous functions on $[0,1]$, determined by linear restrictions on $f$ and possibly its higher-order derivatives. Common examples include bounded, monotone, convex, and concave functions.

As discussed in the introduction, a general approach is to expand $f$ in some basis $\{\phi_{j}\}$ as $f(\cdot)=\sum_{j=1}^{N}\theta_{j}\phi_{j}(\cdot)$ so that the restrictions on $f$ can be posed as linear restrictions on the vector of basis coefficients $\theta\in\mathbb{R}^{N}$, with the parameter space $\mathcal{C}$ for $\theta$ of the form $\mathcal{C}=\{\theta\in\mathbb{R}^{N}\,:\,A\theta\geq b\}$. For example, when $\mathcal{C}_{f}$ corresponds to monotone increasing functions, the set $\mathcal{C}$ is of the form $\{\theta_{1}\leq\theta_{2}\ldots\leq\theta_{N}\}$ under the Bernstein polynomial basis (Curtis and Ghosh,, 2011) and $[0,\infty)^{N}$ under the integrated triangular basis of Maatouk and Bay, (2017). For sake of concreteness, we shall henceforth work with $\mathcal{C}=[0,\infty)^{N}$. Under such a basis representation, the model (3.1) can be expressed as

where $Y=(y_{1},\ldots,y_{n})^{\mathrm{\scriptscriptstyle T}}$ and $\Phi=\{\phi_{j}(x_{i})\}_{ij}$ is an $n\times N$ basis matrix.

The truncated normal prior $\theta\sim\mathcal{N}_{\mathcal{C}}(\mathbf{0},\Omega_{N})$ is conjugate, with the posterior $\theta\mid Y\sim\mathcal{N}_{\mathcal{C}}(\mu_{N},\Sigma_{N})$, with

To motivate their prior choice, Maatouk and Bay, (2017) begin with an unconstrained mean-zero Gaussian process prior on $f$, $f\sim\textsc{gp}(0,K)$, with covariance kernel $K$. Since their basis coefficients correspond to evaluation of the function and its derivatives at the grid points; see Appendix A for details; this induces a multivariate zero-mean Gaussian prior $\mathcal{N}(\mathbf{0},\Sigma_{N})$ on $\theta$ provided the covariance kernel $K$ of the parent Gaussian process is sufficiently smooth. Having obtained this unconstrained Gaussian prior on $\theta$, Maatouk and Bay, (2017) multiply it with the indicator function $\mathbbm{1}_{\mathcal{C}}(\theta)$ of the truncation region to obtain the truncated normal prior.

We are now in a position to connect the posterior bias in Figures 1 and 2 to the mass-shifting phenomenon characterized in the previous section. Since the posterior $\theta\mid Y\sim\mathcal{N}_{\mathcal{C}}(\mu_{N},\Sigma_{N})$, a draw from the posterior can be represented as

Consider an extreme scenario when the true function is entirely flat. In this case, the optimal parameter value $\theta_{0}=\mathbf{0}_{N}$ and under mild assumptions, $\mu_{N}$ is concentrated near the origin with high probability under the true data distribution; see §E of the supplementary material. The mass shifting phenomenon pushes $\theta_{c}$ away from the origin, resulting in the bias. On the other hand, when the true function is strictly monotone as in the left panel of Figure 1, all the entries of $\mu_{N}$ are bounded away from zero, which masks the effect of the shift in $\theta_{c}$.

In strict technical terms, our theory is not directly applicable to $\theta_{c}$ since the scale matrix $\Sigma_{N}$ is a dense matrix in general. However, we show below that $\Sigma_{N}$ is approximately banded under mild conditions. Figure 6 shows image plots of $\Sigma_{N}$ for three choices of $N$ using the basis of Maatouk and Bay, (2017) and sample size $n=500$. In all cases, $\Sigma_{N}$ is seen to have a near-banded structure.

We make this empirical observation concrete below. We first state the assumptions on the basis matrix $\Phi$ and prior covariance matrix $\Omega_{N}$ that allows the construction of a strictly banded matrix approximation.

One example of a basis satisfying Assumption 1 is a B-Spline of fixed order $q$ denoted as $B_{N,q}(x)$ with $N=J+q$ over quasi-uniform knot points of number $J>0$; see, for example, Yoo et al., (2016).

Regarding the prior covariance matrix, we first define a uniform class of symmetric positive definite well-conditioned matrices (Bickel et al., 2008b, ) as

Assumption 2 ensures the covariance matrix is “approximately bandable”, which is common in covariance matrix estimation with thresholding techniques (Bickel et al., 2008a, ; Bickel et al., 2008b, ).

Given above Assumptions, we are now ready to give the approximation result of posterior scale matrix $\Omega^{-1}$ to a banded symmetric positive definite matrix.

Proposition 3 states under mild conditions we can always construct a banded positive definite matrix that approximates $\Sigma_{N}$ in operator norm. Applying the result in Theorem 2 to a truncated normal distribution with the banded approximation of the posterior scale matrix, the marginal density would present a mass-shifting behavior. If we control the band width $K$ such that the approximation is close enough to the posterior scale matrix, the marginal posterior distribution would be expected to behave similarly and shift its probability mass away from the origin and the probability mass over the “corner region” will decrease to zero as the dimension $N$ goes to infinity. This helps explain the bias occurred in the posterior mean over the flat area shown in the Figure 1.

As concluded in previous sections, we view the mass-shifting behavior of the posterior marginals causing the bias in posterior estimation for flat functions. In this section, we will provide empirical evidences that a simple modification to the truncated normal prior can alleviate the issues related to such mass-shifting phenomenon. Among remedies proposed in the literature, Curtis and Ghosh, (2011) proposed independent shrinkage priors on the parameter vector $\{u_{k}\}$ given by a mixture of a point-mass at zero and a univariate normal distribution truncated to the positive real line as an alternative to the truncated normal prior,

Similar mixture priors were also previously used by Neelon and Dunson, (2004) and Dunson, (2005). The mass at zero allows positive prior probability to functions having exactly flat regions. Although possible in principle, introduction of such point-masses while retaining the dependence structure between the coefficients becomes somewhat cumbersome in addition to being computationally burdensome. With such motivation and the additional consideration that in most real scenarios a function is approximately flat in certain regions, we propose a shrinkage procedure as a remedy to replace the coefficients $\theta\in\mathcal{C}$ by $\xi=(\xi_{1},\ldots,\xi_{N})^{\mathrm{\scriptscriptstyle T}}$, where

The parameter $\tau$ provides global shrinkage towards the origin while the $\lambda_{j}$s provide coefficient-specific deviation. We consider default (Carvalho et al.,, 2010) half-Cauchy priors $\mathcal{C}_{+}(0,1)$ on $\tau$ and the $\lambda_{j}$s independently. The $\mathcal{C}_{+}(0,1)$ distribution has a density proportional to $(1+t^{2})^{-1}\mathbbm{1}_{(0,\infty)}(t)$. We continue to use a dependent truncated normal prior $\theta\sim\mathcal{N}_{\mathcal{C}}(\mathbf{0}_{N},\Sigma_{N})$ which in turn induces dependence among the $\xi_{j}$s. Our prior on $\xi$ can thus be considered as a dependent extension of the global-local shrinkage priors (Carvalho et al.,, 2010) widely used in the high-dimensional regression context. Figure 9 in §D.1 of the supplementary material shows prior draws for the first and third components of both $\theta$ and $\xi$, based on which the marginal distribution of the $\xi_{j}$s is clearly seen to place more mass near the origin while retaining heavy tails.

We provide an illustration of the proposed shrinkage procedure in the context of estimating monotone functions as described in (A.1). The procedure can be readily adapted to include various other constraints. Replacing $\theta$ by $\xi$ in (M) in (A.1), we can write (3.1) in vector notation as

Here, $\Psi$ is an $n\times N$ basis matrix with $i^{th}$ row $\Psi_{i}^{\mathrm{\scriptscriptstyle T}}$ where $\Psi_{ij}=\psi_{j-1}(x_{i})$ for $j=1,\ldots,N$ and the basis functions $\psi_{j}$ are as in (A.1). Also, $Y=(y_{1},\ldots,y_{n})^{\mathrm{\scriptscriptstyle T}}$, $\Lambda=\mbox{diag}(\lambda_{1},\ldots,\lambda_{N})$ and $\varepsilon=(\epsilon_{1},\ldots,\epsilon_{n})^{\mathrm{\scriptscriptstyle T}}$.

The model is parametrized by $\xi_{0}\in\mathbb{R}$, $\theta=(\theta_{1},\ldots,\theta_{N})^{\mathrm{\scriptscriptstyle T}}\in\mathcal{C}$, $\lambda=(\lambda_{1},\ldots,\lambda_{N})^{\mathrm{\scriptscriptstyle T}}\in\mathcal{C}$, $\sigma\in\mathbb{R}^{+}$ and $\tau\in\mathbb{R}^{+}$. We place a flat prior $\pi(\xi_{0})\propto 1$ on $\xi_{0}$. We place a truncated normal prior $\mathcal{N}_{\mathcal{C}}(\mathbf{0}_{N},\Sigma_{N})$ on $\theta$ independently of $\xi_{0}$, $\tau$ and $\lambda$ with

and $k(\cdot)$ is the stationary Matérn kernel with smoothness parameter $\nu>0$ and length-scale parameter $\ell>0$. To complete the prior specification, we place improper prior $\pi(\sigma^{2})\propto 1/\sigma^{2}$ on $\sigma^{2}$ and compactly supported priors $\nu\sim\mathcal{U}(0.5,1)$ and $\ell\sim\mathcal{U}(0.1,1)$ on $\nu$ and $\ell$. We develop a data-augmentation Gibbs sampler which combined with the embedding technique of Ray et al., (2019) results in an efficient MCMC algorithm to sample from the joint posterior of $(\xi_{0},\theta,\lambda,\sigma^{2},\tau^{2},\nu,\ell)$; the details are deferred to §D.2 of the supplementary material.

We conduct a small-scale simulation study to illustrate the efficacy of the proposed shrinkage procedure. We consider model (3.1) with true $\sigma=0.5$ and two different choices of the true $f$, namely,

for $x\in[0,1]$. The function $f_{1}$, which is non-decreasing and flat between $0$ and $0.6$, was used as the motivating example in the introduction. The function $f_{2}$ is also approximately flat between $0.7$ and $1$, although it is not strictly non-decreasing in this region, which allows us to evaluate the performance under slight model misspecification.

To showcase the improvement due to the shrinkage, we consider a cascading sequence of priors beginning with only a truncated normal prior and gradually adding more structure to eventually arrive at the proposed shrinkage prior. Specifically, the variants considered are
No shrinkage and fixed hyperparameters: Here, we set $\Lambda=\mathrm{I}_{N}$ and $\tau=1$ in (4.2), and also fix $\nu$ and $\ell$, so that we have a truncated normal prior on the coefficients. This was implemented as part of the motivating examples in the introduction. We fix $\nu=0.75$ and $\ell$ so that the correlation $k(1)$ between the maximum separated points in the covariate domain equals $0.05$.
No shrinkage with hyperparameter updates: The only difference from the previous case is that $\nu$ and $\ell$ are both assigned priors described previously and updated within the MCMC algorithm.
Global shrinkage: We continue with $\Lambda=\mathrm{I}_{N}$ and place a half-Cauchy prior on the global shrinkage parameter $\tau$. The hyperparameters $\nu$ and $\ell$ are updated.
Global-local shrinkage: This is the proposed procedure where the $\lambda_{j}$s are also assigned half-Cauchy priors and the hyperparameters are updated.

We generate $500$ pairs of response and covariates and randomly divide the data into $300$ training samples and $200$ test samples. For all of the variants above, we set the number of knots $N=150$. We provide plots of the function fit along with pointwise $95\%$ credible intervals in Figures 7 and 8 respectively, and also report the mean squared prediction error (mspe) at the bottom of the sub-plots. As expected, only using the truncated normal prior leads to a large bias in the flat region. Adding some global structure to the truncated normal prior, for instance, updating the gp hyperparamaters and adding a global shrinkage term improves estimation around the flat region, which however still lacks the flexibility to transition from the flat region to the strictly increasing region. The global-local shrinkage performs the best, both visually and also in terms of mspe.

Additionally, the shrinkage procedure performed at least as good as bsar, a very recent state-of-the-art method, developed by Lenk and Choi, (2017), and implemented in the R package bsamGP. For the out-of-sample prediction performance of bsar, refer to Figure 10 in §D.3 of the supplementary material, based on which it is clear that the performance of global-local shrinkage procedure is comparable with that of bsar. It is important to point out that bsar is also a shrinkage based method that allows for exact zeros in the coefficients in a transformed Gaussian process prior through a spike and slab specification.

A seemingly natural way to define a prior distribution on a constrained parameter space is to consider the restriction of a standard unrestricted prior to the constrained space. The conjugacy properties of the unrestricted prior typically carry over to the restricted case, facilitating computation. Moreover, reference priors on constrained parameters are typically the unconstrained reference prior multiplied by the indicator of the constrained parameter space (Sun and Berger,, 1998). Despite these various attractive properties, the findings of this article pose a caveat towards routine truncation of priors in moderate to high-dimensional parameter spaces, which might lead to biased inference. This issue gets increasingly severe with increasing dimension due to the concentration of measure phenomenon (Talagrand,, 1995; Boucheron et al.,, 2013), which forces the prior to increasingly concentrate away from statistically relevant portions of the parameter space. A somewhat related issue with certain high-dimensional shrinkage priors has been noted in Bhattacharya et al., (2016). Overall, our results suggest a careful study of the geometry of truncated priors as a useful practice. Understanding the cause of the biased behavior also suggests natural shrinkage procedures that can guard against such unintended consequences. We note that post-processing approaches based on projection (Lin and Dunson,, 2014) and constraint relaxation (Duan et al.,, 2020) do not suffer from this unintended bias. The same is also true for the recently proposed monotone bart (Bayesian Additive Regression Trees) (Chipman et al.,, 2010) method.

It would be interesting to explore the presence of similar issues arising from truncations beyond the constrained regression setting. Possible examples include correlation matrix estimation and simultaneous quantile regression. Priors on correlation matrices are often prescribed in terms of constrained priors on covariance matrices, and truncated normal priors are used to maintain ordering between quantile functions corresponding to different quantiles, and this might leave the door open for unintended bias to creep in.