Truncated multivariate normal

[papamarkou]

I am interested in implementing the quantile() method for MvNormal, as I need to be able to compute the quantiles of MvNormal (so that I can ultimately implement a truncated multivariate normal distribution similar in spirit to the truncated univariate normal).

I have tried to find out what is the most reasonable and efficient way to do this, and it seems there is a closed form for it here. I will try to look into the R code of qmvnorm as well, which uses an algorithm by Genz and Bretz, see this doc page.

I will try to find out more about this matter. If any of you any suggestions about the choice of algorithm or implementation, I am open to any new thoughts.

[papamarkou]

I substantially changed the title of the issue, because after inspecting more carefully my goal of defining a truncated multivariate normal, I realized that the quantile computation is not necessarily so central to the issue.

Furthermore, after looking at the code in Distributions/src/truncate.jl, it appears that the Truncated type is a sub-type of UnivariateDistribution, which doesn't make sense :) Why such restriction? Truncated distributions can be as well multivariate, see for example this nice exposition of truncated MvNormals.

[johnmyleswhite]

I'd be cautious about the R code for qmvnorm since that package is GPL. If you've already read that code, you probably can't write an equivalent of that method for Distributions.jl -- it will have to go into a GPL package instead.

We could create an abstract AbstractTruncated type and/or rename the current type, which is specialized to the univariate setting to allow us to implement it using a simple concrete type. It's probably just a reflection of my lack of knowledge, but I don't see how to handle generic truncation in multiple dimensions, so I don't know how to make a generic type for truncation outside of the univariate case. If you figure that out, there's probably no need to add any more layers of abstract types.

[papamarkou]

I am not sure either what would be the correct definition of a mutlivariate truncated distribution. I am thinking out loud - would the following definition of MvTruncated be legit?

immutable MvTruncated{D<:MultivariateDistribution, S<:ValueSupport} <: MultivariateDistribution{S}
    untruncated::D      # the original distribution (untruncated)
    lower::Vector{Float64}      # lower bound
    upper::Vector{Float64}      # upper bound
    lcdf::Float64       # cdf of lower bound
    ucdf::Float64       # cdf of upper bound

    tp::Float64         # the probability of the truncated part, i.e. ucdf - lcdf
    logtp::Float64      # log(tp), i.e. log(ucdf - lcdf)
end

All I changed in comparison to the univariate case is that now the lower and upper fields are vectors.

[johnmyleswhite]

It seems like you're going to end up recreating the kinds of constraints that people use in optimization: here you've got box constraints. Seems reasonable, but I have no idea what people are using truncated multivariate Gaussians for.

[papamarkou]

To give one simple example, let's say that I have a pair of (correlated) parameters (u, v) I want to sample from using MCMC, and let's further say that u and v are constrained. If I want to use Metropolis-Hastings to sample (u, v), I can define a truncated bivariate normal as the proposal at each iteration of the sampler, centered at the state of the current iteration (additionally, I would need to correct the acceptance ratio by the normalising constant that relates to the cdf of the underlying multivariate normal).

A simpler approach would be to carry on sampling from the multivariate normal till I produce (u, v) that satisfy the constraint; this would be the same as sampling from the truncated multivariate normal, only it may be less efficient (in this case the normalising constant correction in the ratio would still be needed).

What I don't know yet is how easy it is to compute the normalising constant (cdf) of multivariate normal, i.e. F(upper)-F(lower), where now lower and upper are vectors.

[fiatflux]

I just ran into a perfectly justifiable scientific case for multivariate normal distributions truncated to the cone with positive entries. High dimensionality is the main reason a collection of univariates would be unsatisfactory, although I haven't verified the present implementation gains us a whole lot to that end.

[matbesancon]

This seems rather complex and doable outside of Distributions, closing?

[filchristou]

Since Product was introduced in 2018, it's worth mentioning that today this can be (partially) done with Product([truncated(Normal(μ,σ), low, high) for (μ,σ) in zip(μs,σs)]) but still only for independent distributions.