<p>I have reimplemented an algorithm which does not depend on MCMC but creates <strong>independent and identically distributed (iid) samples</strong> from the truncated multivariate normal distribution. Having iid samples can be very useful! I used to also use <a href="https://emcee.readthedocs.io/en/stable/" rel="nofollow noreferrer">emcee</a> as described in the answer by Warrick, but for convergence the number of samples needed exploded in higher dimensions, making it impractical for my use case.</p>
<p>The algorithm was introduced by <a href="https://arxiv.org/pdf/1603.04166.pdf" rel="nofollow noreferrer">Botev (2016)</a> and uses an accept-reject algorithm based on minimax exponential tilting. It was <a href="https://de.mathworks.com/matlabcentral/fileexchange/53792-truncated-multivariate-normal-generator" rel="nofollow noreferrer">originally implemented in MATLAB</a> but reimplementing it for Python increased the performance significantly compared to running it using the MATLAB engine in Python. It also works well and is fast at higher dimensions.</p>
<p>The code is available at: <a href="https://github.com/brunzema/truncated-mvn-sampler" rel="nofollow noreferrer">https://github.com/brunzema/truncated-mvn-sampler</a>.</p>
<h2>An Example:</h2>
<pre><code>d = 10  # dimensions

# random mu and cov
mu = np.random.rand(d)
cov = 0.5 - np.random.rand(d ** 2).reshape((d, d))
cov = np.triu(cov)
cov += cov.T - np.diag(cov.diagonal())
cov = np.dot(cov, cov)

# constraints
lb = np.zeros_like(mu) - 1
ub = np.ones_like(mu) * np.inf

# create truncated normal and sample from it
n_samples = 100000
tmvn = TruncatedMVN(mu, cov, lb, ub)
samples = tmvn.sample(n_samples)
</code></pre>
<p>Plotting the first dimension results in:</p>
<p><a href="https://i.stack.imgur.com/WxWxq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/WxWxq.png" alt="First dimension of truncated MVN" /></a></p>
<hr />
<h2>Reference:</h2>
<p>Botev, Z. I., (2016), The normal law under linear restrictions: simulation and estimation via minimax tilting, Journal of the Royal Statistical Society Series B, 79, issue 1, p. 125-148</p>
Truncated multivariate normal in SciPy?
<p>I'm trying to automate a process that at some point needs to draw samples from a truncated multivariate normal.  That is, it's a normal multivariate normal distribution (i.e. Gaussian) but the variables are constrained to a cuboid.  My given inputs are the mean and covariance of the full multivariate normal but I need samples in my box. </p>

<p>Up to now, I'd just been rejecting samples outside the box and resampling as necessary, but I'm starting to find that my process sometimes gives me (a) large covariances and (b) means that are close to the edges.  These two events conspire against the speed of my system.</p>

<p>So what I'd like to do is sample the distribution correctly in the first place.  Googling led only to <a href="https://groups.google.com/forum/#!topic/scipy-user/QohTys9U06M" rel="noreferrer">this discussion</a> or the <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.truncnorm.html#scipy.stats.truncnorm" rel="noreferrer"><code>truncnorm</code> distribution</a> in <code>scipy.stats</code>.  The former is inconclusive and the latter seems to be for one variable.  Is there any native multivariate truncated normal?  And is it going to be any better than rejecting samples, or should I do something smarter?</p>

<p>I'm going to start working on my own solution, which would be to rotate the untruncated Gaussian to it's principal axes (with an SVD decomposition or something), use a product of truncated Gaussians to sample the distribution, then rotate that sample back, and reject/resample as necessary.  If the truncated sampling is more efficient, I think this should sample the desired distribution faster.</p>