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