将 Julia 中的正态分布与 Distributions.jl 结合起来的正确方法

Question

我有两个这样的多元正态分布：

using Distributions, LinearAlgebra
g1 = MvNormal([1,2], [2 1; 1 2])
A = [3 1; 1 3]
B = [[A [0;0]]; transpose([0,0,1])]
g2 = MvNormal([1,2,3], B)

我想将它们组合成两个变量串联的新分布，假设它们是独立的。函数 Product_distribution 似乎应该可以解决问题：

g3 = product_distribution(g1, g2)

但这会导致错误：

ERROR: all distributions must be of the same size

我真的不明白，这让我觉得这个函数的目的不是我想象的那样，但我找不到任何其他更合适的。

需要明确的是，所需的输出应等于：

m3 = vcat(mean(g1), mean(g2))
s3 = hvcat( (2,2), cov(g1), zeros(2,3), zeros(3,2), cov(g2))
g3 = MvNormal(m3, s3)

（虽然也许稀疏矩阵或其他优化的对角块矩阵类型会更合适，但在这种情况下我真的不在乎。）

Answer 1

找不到明确的答案，但这个问题之前已经出现过。到目前为止我能建议的最好的：

using BlockDiagonals

function concat(d1::Union{MvNormal, Normal}, d2::Union{MvNormal, Normal})
    c1 = d1 isa MvNormal ? cov(d1) : [var(d1);;]
    c2 = d2 isa MvNormal ? cov(d2) : [var(d2);;]
    MvNormal(vcat(mean(d1),mean(d2)), BlockDiagonal([c1, c2]))
end

给予：

julia> g1 = MvNormal([1,2], [2 1; 1 2]);

julia> B = BlockDiagonal([A,[1;;]]);

julia> g2 = MvNormal([1,2,3], B);

julia> g12 = concat(g1,g2)
MvNormal{Float64, PDMats.PDMat{Float64, BlockDiagonal{Float64, Matrix{Float64}}}, Vector{Float64}}(
dim: 5
μ: [1.0, 2.0, 1.0, 2.0, 3.0]
Σ: [2.0 1.0 … 0.0 0.0; 1.0 2.0 … 0.0 0.0; … ; 0.0 0.0 … 3.0 0.0; 0.0 0.0 … 0.0 1.0]
)

julia> rand(g12)
5-element Vector{Float64}:
  1.0362522596298223
  2.0469956784329866
 -0.21925320262748982
 -1.2334419613775114
  3.2146164519549814