{'fraction_non_alphanumeric': 0.04511182108626198, 'fraction_numerical': 0.017678381256656016, 'mean_word_length': 4.282628262826282, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Soft Multivariate Truncated Normal Distribution

We propose a new distribution, called the soft tMVN distribution, which provides a smooth approximation to the truncated multivariate normal (tMVN) distribution with linear constraints. An efficient blocked Gibbs sampler is developed to sample from the soft tMVN distribution in high dimensions. We provide theoretical support to the approximation capability of the soft tMVN and provide further empirical evidence thereof. The soft tMVN distribution can be used to approximate simulations from a multivariate truncated normal distribution with linear constraints, or itself as a prior in shape-constrained problems.


Introduction
The truncated multivariate normal (tMVN) distribution is routinely used as a prior distribution on model parameters in Bayesian shape-constrained regression. Structural constraints, such as monotonicity and/or convexity, are commonly induced by expanding the function in an appropriate basis where the constraints can be induced by imposing linear constraints on the coefficients; some examples of such a basis include piecewise linear functions [10], splines [6], Bernstein polynomials [34], and compactly supported basis functions [23]. Under a Gaussian or scale-mixture of Gaussian error distribution, the conditional posterior of the basis coefficients once again turns out to be truncated normal with linear constraints [23].
The problem of sampling from a tMVN distribution with linear constraints is also frequently encountered as a component of a larger Markov chain Monte Carlo (MCMC) algorithm to sample from the full conditional distribution of a constrained parameter vector. As a running example revisited on multiple occasions in this article, consider binary variables y i = 1(z i > 0), with z = (z 1 , . . . , z n ) T a vector of latent Gaussian thresholds [1] and w ∈ R q a vector of parameters/latent variables so that the joint distribution of θ = (z, w) follows a N (µ, Σ) distribution. It then immediately follows that the (conditional) posterior of θ | y, µ, Σ follows a N (µ, Σ) distribution truncated to ⊗ n i=1 C i ⊗ R q , with C i = (0, ∞) or (−∞, 0) depending on whether y i = 1 or 0. Such latent Gaussian threshold models are ubiquitous in the analysis of binary and nominal data; examples include probit regression and its multivariate extensions [1,17,7,26], multinomial probit models [24,36,18], tobit models [33,28], and binary Gaussian process (GP) classification models [16] among others.
In this article, we propose a new family of distributions called the soft tMVN distribution which replaces the hard constraints in a tMVN distribution with a smoothed or "soft" version using a logistic sigmoid function. The soft tMVN distribution admits a smooth log-concave density on the d-dimensional Euclidean space. Although the soft tMVN distribution is supported on the entire d-dimensional space, it can be made to increasingly concentrate most of its mass on a polyhedron determined by multiple linear inequality constraints, by tweaking a parameter. In fact, we show that the soft tMVN distribution approximates the corresponding tMVN distribution in total variation distance.
Recognizing the soft tMVN distribution as the posterior distribution in a pseudo-logistic regression model, we develop an efficient blocked Gibbs sampler combining the Polya-Gamma data augmentation of Polson, Scott & Windle (2013) [30] along with a structured multivariate normal sampler from Bhattacharya, Chakraborty, and Malllick (2016) [3]. In contrast, existing Gibbs samplers for a tMVN distribution sample the coordinates one-at-a-time from their respective full conditional univariate truncated normal distributions [14,20,9,32]. The algorithm of Geweke is implemented in the R package tmvtnorm [35]. While the Gibbs sampling procedure is entirely automated, it is well-recognized in a broader context that such one-at-a-time updates can lead to slow mixing, especially if the variables are highly correlated. We have additionally observed numerical instabilities in the R implementation for unconstrained dimensions exceeding 400. While exact Hamiltonian Markov chain (HMC) algorithms to sample from tMVN [27] are also popular, such algorithms are not suitable to sample from the soft tMVN, and leaf-frog steps with careful tuning are necessary to obtain good mixing. There also exists accept-reject algorithms for the tMVN distribution that create exact samples from the distribution [4]. The algorithm of Botev is implemented in the R package TruncatedNormal [5]. While exact samples are possible, when the acceptance probability becomes small, either the algorithm slows tremendously or approximate samples are produced. We typically saw small acceptance probabilities in the R implementation when the constrained dimension exceeded 200. With such motivation, we propose to replace a tMVN distribution with its softened version inside a larger MCMC algorithm and use our sampling strategy for the soft tMVN distribution. In recent years, there has been several instances of such approximate MCMC (aMCMC) [19] algorithms where the exact transition kernel of a Markov chain is replaced by an approximation thereof for computational ease.
The soft tMVN distribution can also be used as a prior distribution in Bayesian shape-constrained regression problems as an alternative to the usual tMVN prior. Like the tMVN distribution, the soft tMVN distribution is conditionally conjugate for the mean in a Gaussian likelihood. The soft tMVN can be viewed as a shrinkage prior which encourages shrinkage towards a linearly constrained region rather than being supported on the region. There is an interesting parallel between the soft tMVN distribution and global-local shrinkage priors used in sparse regression problems. The global-local priors replace the point mass (at zero) of the more traditional discrete mixture priors and rather encourage shrinkage towards the origin, with the motivation that a subset of the regression coefficients may have a small but non-negligible effect. Similarly, the soft tMVN prior favors the shape constraints while allowing for small departures.
The rest of the article is organized as follows. In Section 2, we introduce the soft tMVN distribution as an approximation to the tMVN distribution and discuss its properties. In Section 3, we discuss various strategies to sample from a soft tMVN distribution, including a scalable Gibbs sampler suitable for high-dimensional situations. Section 4 contains a number of simulation examples to illustrate the efficacy of the proposed sampler as well as the approximation capability of the soft tMVN distribution. We conclude with a discussion in Section 5.

The soft tMVN distribution
Consider a tMVN distribution where µ ∈ R d , Σ is a d × d positive definite matrix, and C is described by r ≤ d linear constraints, where s i ∈ {1, −1} denotes the sign of the ith inequality, and a i ∈ R d . Without loss of generality, we assume the first r coordinates to be constrained; this is mainly for notational convenience and can always be achieved by reordering the variables, if necessary. We also assume throughout that C has positive R d -Lebesgue measure, so that the density γ in (1) is non-singular on R d . In the special case where a i = e i , the ith unit vector in R d (with 1 at the ith coordinate and 0 elsewhere), the constraint set C reduces to the form ⊗ r i=1 C i ⊗ R q mentioned in the introduction. While this is an important motivating example, our approach works more generally for the type of constraints in the above display.
Write, using the convention 0 0 = 1, Our main idea is to replace the indicator functions above with a smoothed or "soft" approximation.
A rich class of approximations to the indicator function 1 (0,∞) (·) is provided by sigmoid functions, which are non-negative, monotone increasing, differentiable, and satisfy lim x→∞ σ(x) = 1 and lim x→−∞ σ(x) = 0. The cumulative distribution function of any absolutely continuous distribution on R which is symmetric about zero can be potentially used as a sigmoid function. Here, for reasons to be apparent shortly, we choose to use the logistic sigmoid function σ(x) = 1/(1 + e −x ), which is the cdf of the logistic distribution. Specifically, define, for η > 0, to be a scaled version of σ(·). The parameter η controls the quality of the approximation, with larger values of η providing increasingly better approximations to 1 (0,∞) (·). In fact, it is straightforward to see that It is also immediate that (1 − σ η (·)) is an approximation to 1 (−∞,0) (·) with the same approximation error.
We are now ready to describe our approximation scheme. Fixing some large η and replacing the indicators by their respective sigmoidal approximations in (1), we obtain the approximation γ η to γ as We refer to γ η as a soft tMVN distribution and generically denote it by N s C (µ, Σ). It is immediate to note that γ η is a smooth (infinitely differentiable) density supported on R d . Further, a simple calculation shows that i.e., the Hessian matrix of the negative log density is positive definite. This implies that γ η is a logconcave density, which, in particular means γ η is unimodal. We collect these various observations about γ η in Proposition 2.1.
Proposition 2.1. Let γ and γ η be respectively defined as in (1) and (4). Then, γ η is an infinitely differentiable, unimodal, log-concave density on R d . Further, A proof is provided in the Appendix. The last part of Proposition 2.1 formalizes the intuition that γ η approximates γ for large η by showing that the L 1 distance between γ η and γ converges to 0 as η → ∞. An inspection of the proof for the L 1 approximation will reveal that we haven't used any particular feature of the logistic function and the argument can be extended to other sigmoid functions.
The L 1 approximation result implies that although γ η has a non-zero density at all points in R d , the effective support is the region C for large values of η, and a random draw from γ η will fall inside C with overwhelmingly large probability. To obtain a more quantitative feel for how the approximation gets better with increasing η, we set γ to be a standard bivariate normal distribution truncated to the first orthant, Figure 1 shows contour plots of γ (last column) along with those for γ η for various values of η, with η increasing from left to right. Each row corresponds to a different value of ρ. It is evident that the approximation quickly improves at η increases, and stabilizes around η = 100. We later show in simulations involving substantial higher dimensions that γ η with η = 100 continues to provide a reasonable approximation to the corresponding tMVN distribution γ. The accurate approximation of the soft tMVN has two important consequences in our opinion. First, for any of the examples discussed in the introduction which require a sample from a tMVN within an MCMC algorithm, a sample from a tMVN can be replaced with a sample from the corresponding soft tMVN distribution; we discuss efficient strategies to sample the tMVN distribution in the next section. Second, the soft tMVN distribution can itself be used as a prior distribution for constrained parameters. As a prior, the soft tMVN replaces the hard constraints imposed by the tMVN with soft constraints, encouraging shrinkage towards the constrained region C. Indeed, the soft tMVN distribution can be considered a global shrinkage prior [29] which shrinks vectors towards a pre-specified constrained region.
The tMVN prior is conditionally conjugate for a Gaussian likelihood and the soft tMVN prior naturally inherits this conditional conjugacy. Suppose Y | θ, σ 2 ∼ N (Φθ, σ 2 I n ) and θ ∼ N s C (µ, Σ) is assigned a soft tMVN prior. Then, The conditional conjugacy allows one to fit a conditionally Gaussian model with a soft tMVN prior using standard Gibbs sampling algorithms, provided one can efficiently sample from a soft tMVN distribution. We discuss this in the next section.
3 Sampling from the soft tMVN distribution

Gibbs sampler in high-dimensions
In this subsection, we propose a scalable data-augmentation blocked-Gibbs sampler to sample from a soft tMVN distribution. The proposed Gibbs sampler updates the entire θ vector in a block, unlike one-at-a-time updates for Gibbs samplers for tMVNs. Apart from log-concavity, the other nice feature behind our choice of the logistic sigmoid function is that γ η can be recognized as the posterior distribution of a vector of regression parameters in a logistic regression model. To see this, consider the setup of a logistic regression model with binary response t i ∈ {0, 1} and vector of predictors W i ∈ R d for i = 1, . . . r, Assuming a N (µ, Σ) prior on the vector of regression coefficients θ, the posterior distribution of θ | t, W, µ, Σ is given by If we now set t i = 1(s i = 1) and W i = η a i , then the above density is identical to γ η . The number of constraints r plays the role of the sample size, and the ambient dimension d ≥ r indicates the number of the regression parameters in this pseudo-logistic model. Thus, sampling from γ η is equivalent to sampling from the conditional posterior of regression parameters in a high-dimensional logistic regression model, which can be conveniently carried out using the Polya-Gamma data augmentation scheme of Polson, Scott & Windle (2013) [30]. The Polya-Gamma scheme introduces r auxiliary variables ω 1 , . . . , ω r and performs Gibbs sampling by alternatively sampling from ω | θ, t and θ | ω, t as follows: where W ∈ R r×d with ith row W i , t = (t 1 , . . . , t r ) , κ = (t − 1/2), and Ω = diag(ω 1 , . . . , ω r ).
In (i), PG denotes a Polya-Gamma distribution which can be sampled using the Bayeslogit package in R [30]. Note that the entire θ vector is sampled in a block in step (ii). The worstcase complexity of sampling from the multivariate Gaussian distribution in (6) is O(d 3 ). However, exploiting the structure of µ ω and Σ ω , a sample from N (µ ω , Σ ω ) can be obtained with significantly less cost using a recent algorithm in Bhattacharya et al. (2016) [3] provided d r and a N (0, Σ) variate can be cheaply sampled.
Define Φ = Ω 1/2 W and α = Ω −1/2 κ. Then, a sample from (ii) is obtained by first samplinḡ and setting First, by the Sherman-Woodbury-Morrison formula, which only requires solving a r × r system. Samplingθ in (7)   where L is a q × q diagonal matrix and H is an N × q (possibly dense) matrix. A sample u from N (0, Σ) is then obtained by (i) Sample z ∼ N (0, I N ) and u 2 ∼ N (0, L) independently.
(ii) Set u 1 = Hu 2 + z and u = (u 1 , u 2 ) . Since u is a linear transformation of (z, u 2 ) which is jointly Gaussian, u also has a joint Gaussian distribution. Calculating the covariance matrix of u then immediately shows that u ∼ N (0, Σ). Since L is diagonal, u 2 can be sampled in O(q) steps, and the matrix multiplication costs O(N q 2 ), so that the overall cost is O(N q 2 ).

Other strategies
In moderate dimensions, it is possible to use a Metropolis (Gaussian) random walk and its various extensions to sample from a soft tMVN distribution. In particular, given that the soft tMVN distribution can be recognized as the posterior distribution in a model with a Gaussian prior, elliptical slice sampling [25] is a viable option. There is substantial literature on sampling from log-concave distributions using variants of the Metropolis algorithm with strong theoretical guarantees [13,12,21,22,2]. More recently, Dalalyan (2017) [8] and Durmus & Moulines (2016) [11] provided non-asymptotic bounds on the rate of convergence of unadjusted Langevin Monte Carlo (LMC) algorithms for log-concave target densities. Assuming the target density is proportional to e −f (θ) for some convex function f , the successive iterates of a first-order LMC algorithm takes the form where the {ξ k }s are independent N (0, I) variates and h > 0 is a step-size parameter. Clearly, {θ k } k=0,1,... forms a discrete-time Markov chain and the results in Dalalyan (2017) [8] and Durmus & Moulines (2016) [11] characterize the rate at which the distribution of θ k converges to the target density in total variation distance. Aside from the non-asymptotic bounds, another key message from their results is that the typical Metropolis adjustment as in Metropolis adjusted Langevin (MALA) [?] is not required for log-concave targets. Dalalyan (2017) [8] also provide a second-order version of the LMC algorithm called LMCO which can incorporate the Hessian ∇ 2 f . Since both ∇(− log γ η ) and ∇ 2 (− log γ η ) are analytically tractable, it is possible to use both the LMC and LMCO algorithms to sample from γ η . Other than MCMC, another possible strategy to sample from γ η is to use a multivariate generalization of the adaptive rejection sampling (ARS) [15].

Simulations
In this section, we conduct a number of simulations to empirically illustrate that the soft tMVN distribution continues to provide an accurate approximation to the tMVN distribution in highdimensional situations. These simulations also demonstrate the scalability of the proposed Gibbs sampler.
To begin with, we first justify our continued use of η = 100 in higher dimensions. In Figure 1, we had provided the contour plots of a bivariate tMVN distribution and its soft tMVN approximation with η = 100. As an obvious extension, we now consider the bivariate marginal of (θ 1 , θ 2 ), where θ ∈ R 50 is drawn from a multivariate normal distribution with mean µ = 0 and with a compound symmetry covariance structure, Σ = (1 − ρ)I 50 + ρ1 50 1 T 50 , truncated to the positive orthant. We consider two choices of ρ, namely ρ = 0.25 and 0.75, and provide the contour plots for N C (µ, Σ) and N s C (µ, Σ) in the top and bottom panels of Figure 2 respectively. The contour plots were drawn by collecting 150, 000 samples from the N C (µ, Σ) and N s C (µ, Σ) distributions, and then retaining the first two coordinates in each case to obtain samples from the bivariate marginal. Specifically, we used the rejection sampler of Botev (2017) [4] implemented in the R package TruncatedNormal [5] to draw samples from a tMVN distribution and used our data augmentation Gibbs sampler to sample from the soft tMVN distribution. The figure shows that η = 100 remains a reasonable choice in higher dimensions, and we henceforth fix η = 100 throughout.
Next, we provide some numerical summaries in two different settings. Due to the inherent difficulty of comparing two high-dimensional distributions, we will compare the marginal densities. Specifically, given densities f and g on R d with finite mean, we consider two different measures to compare them. The first one uses the 1 st Wasserstein (W 1 ) distance between two distributions, W 1 (f, g). It is an average W 1 distance between the marginals, where f i denotes the ith marginal density of f . Our second measure is an average squared L 2 distance between the mean vectors for the two densities, with µ f = R d xf (x)dx. We compute D and ξ between γ and γ η for two different covariance structures in Σ. Due to the lack of analytic expressions for the marginals for non-diagonal Σ, we resort to simulations to approximate D and ξ. The highest dimension d used in our simulations is d = 600; while our sampler can be scaled beyond this, the rejection sampler starts producing warning messages due to incurring small acceptance probabilities.
is a binary response at locations s i for i = 1, . . . , n. This is modeled as s j ), and K is the Matern kernel. Letting Z = [Z(s 1 ), . . . , Z(s n )] T , the conditional distribution of Z | Y follows the above N n (0, For the simulation, set n = {100, 200}. Let s i = i for i = 1, . . . , n. We randomly sample 1 from {10, . . . , n/2} and 2 from {n/2 + 1, . . . , n − 10} and let Y 1 , . . . , Y 1 = 1, Y 1 +1 , . . . , Y 2 = 0, and Y 2 +1 , . . . , Y n = 1. We set the smoothness parameter for the Matern kernel at 3/5 and the scale parameter at 1. We then proceed to draw 5000 samples from the tMVN, N n (0, Σ)1 C (θ), using Botev's rejection sampler and 5000 samples from the soft tMVN, N s n (0, Σ)1 C (θ), using our Gibbs sampler. The 5000 samples were collected for our method after discarding 1000 initial samples as burn-in and collecting every 100th sample to thin the chain. Figures 3 and 4 show the marginal density plots of 8 coordinates of θ based on the 5000 samples for the two values of n respectively. The tMVN distribution is shown in blue while the soft tMVN is in pink. It is evident that for both values of n, the marginal densities are visually indistinguishable. To obtain an overall summary measure, Figure 5 shows the histogram of ξ, defined in equation (11), (left panel) and D, defined in (10), over 50 independent simulations. Both the histograms are tightly centered near the origin, which again suggests the closeness of the tMVN and soft tMVN distributions. As a quick comparison, the value of D between N(0, Σ) and N(0.005, Σ) for the current Σ is about 0.03 for both values of n.
For this simulation, we sample x i iid ∼ N (0, I P ) and λ j ∼ U [1/15, 1/5], and then set Σ to the above form. Draw β ∼ N (0, Λ) and Z ∼ N (Xβ, I n ). Then if Z i ≥ 0, set Y i = 1 and if Z i < 0, set Y i = 0. For both (N, P ) = {(100, 400), (200, 400)}, we then proceed to draw 5000 samples from the tMVN, N n (0, Σ)1 C (θ), using Botev's rejection sampler and 5000 samples from the soft tMVN, N s n (0, Σ)1 C (θ), using our Gibbs sampler. The 5000 samples were collected for our method after discarding 1000 initial samples as burn-in and collecting every 100th sample to thin the chain. Figures 6 and 7 show the marginal density plots of 8 coordinates of θ based on the 5000 samples for the two combinations respectively; as before, the tMVN distribution is shown in blue while the soft tMVN is in pink. We once again see that for both combinations, the marginal densities overlap well. To obtain an overall summary measure, Figure 8 shows the histogram of ξ, defined in equation (11), (left panel) and D, defined in (10), over 50 independent simulations. We see that the histogram of ξ and D shifts to the right for n = 200 than for n = 100. This shift is expected as the size of the matrix X grows; as a point of comparison, in Figure 9, we plot the histogram of D between N (0, Σ) and N (0.005, Σ) for the present choice of Σ and see a similar shift.

Discussion
In this paper, we have presented the soft tMVN distribution, which provides a smooth approximation to the tMVN distribution with linear constraints. Our theoretical and empirical results suggest that the soft tMVN distribution offers a good approximation to the tMVN distribution in high dimensional situations. Future work will focus on using the soft tMVN distribution as a prior in various shape-constrained models, with posterior computation aided by the Gibbs sampler developed herein. where recall σ η (·) is the sigmoidal function. Bound where we have used triangle inequality and the fact that g η (θ)dθ = M η . Now, we have |M − M η | = (g − g η ) ≤ |g(θ) − g η (θ)|dθ.
Now, using the inequality (2) and the fact that for numbers u i , v i ∈ [0, 1], we have that |g(θ) − g η (θ)|dθ where a b means a ≤ Cb for some positive constant C, and the expectation E is under a N (µ, Σ) distribution. By monotone convergence theorem, the right hand side of the above display converges to 0 as η → ∞.