Hi! Is any package to get Multivariate (normal) Conditional Distribution?

Something like described here .

maybe other people jump in that know more, but this seems related:

opened 10:14PM - 22 Sep 15 UTC

new-func
scope

I'm interested writing a set of functions that condition on a subset of variable… s, and return a new (conditional) distribution. This is _roughly_ how I'd imagine doing it for a multivariate normal distribution:
``` julia
# Given a multivariate distribution p(x),
# returns new distribution p(x | xᵢ=aᵢ), where
# "a" is a vector of values to condition on, and
# "d2" is a vector of indices providing the
# dimensions (subscript i) in increasing order.
function condition(
joint::MultivariateDistribution,
a::Vector{Float64},
d2::Vector{Int}
)
if length(d2) != length(a)
error("list of conditioned values and dimension indices bet the same length")
end
max_dim = length(joint.μ)
if minimum(d2) < 1 || maximum(d2) > max_dim
error("dimension indices must be between 1 and dimension of joint distribution")
end
# d1 = remaining dimensions
# d2 = removed/conditioned dimensions
d2 = sort(d2)
d1 = setdiff(1:max_dim,d2)
# covariance matrix blocks
Σ = full(joint.Σ)
Σ11 = Σ[d1,d1]
Σ12 = Σ[d1,d2]
Σ22inv = inv(Σ[d2,d2])
μ_new = joint.μ[d1] + (Σ12 * Σ22inv * (a - joint.μ[d2]))
Σ_new = Σ11 - (Σ12 * Σ22inv * Σ12')
return MvNormal(μ_new,Σ_new)
end
```
**Questions:**
1) I know this is really inefficient, particularly the step `Σ = full(joint.Σ)` to convert a `PDmat` to a normal matrix. Any suggestions on how to do this better -- i.e. using the (already computed) Cholesky factorization? I admit I haven't given this much thought at all yet.
2) Would a set of functions like this be useful to add to `Distributions.jl` or does it not fit within the scope of this package? Where (if anywhere at all) would this live in JuliaStats?
3) Is the above approach completely off base, or am I on the right track at least?

in particular, also the last answer therein and reference to the `conditional(P, A, B, xB)`

method from GaussianDistributions.jl which should do what you want.

essentially you could “just” calculate \mu_{1|2} and \Sigma_{1|2} yourself as described in your link and then sample from `MvNormal`

(\mu_{1|2} , \Sigma_{1|2} ) via Distributions.jl. Advantage of the above references may be a good numerical implementation for computing \mu_{1|2} and \Sigma_{1|2} .

1 Like