Singular Exception with LKJCholesky

Specific Domains Probabilistic programming
21

Alrighty! It took a bit longer than I originally thought because I realized we should clean up some stuff along the way, but here is the PR: Remove usage of `DynamicPPL.tonamedtuple` by torfjelde · Pull Request #2071 · TuringLang/Turing.jl · GitHub

Once that’s through, LKJCholeksy (and others) should also work nicely with MCMChcains, etc.

1 Like
22

I made 2-dimentional State stace model . But the posterior of the coefficient matrix deviates from what I had inputted. Would you have any insight?

#seed
Random.seed!(1234)
T = 50

B  = [0.9 0.1; 0.0 0.9]
Σ_X = [0.2 0.1; 0.1 0.2]
Σ_Y = [0.1 0.05; 0.05 0.1]
X = zeros(T, 2)
X[1, :] = [0, 0]
Y = zeros(T, 2)
for t in 2:T
    X[t, :] = B * X[t-1, :] + rand(MvNormal(Σ_X))
    Y[t, :] = X[t, :] + rand(MvNormal(Σ_Y))
end

@model function SSM_2D(y, ::Type{TV}=Matrix{Float64}) where {TV}
    T = size(y, 1)
    k = size(y, 2)
    #prior
    ###----
    F_x ~ LKJCholesky(k, 1.0)
    #σ_x : 2 ×1 vector
    σ_x ~ filldist(truncated(Cauchy(0., 5.); lower = 0), k)
    Σ_xL = Diagonal(σ_x .+ 1e-3) * F_x.L
    Σ_x = PDMat(Cholesky(Σ_xL + eps() * I))

    F_y ~ LKJCholesky(k, 1.0)
    σ_y ~ filldist(truncated(Cauchy(0., 5.); lower = 0), k)
    Σ_yL = Diagonal(σ_y .+ 1e-3) * F_y.L
    Σ_y = PDMat(Cholesky(Σ_yL + eps() * I))
    ###---- 
    B = TV(undef, k, k)
    for i in 1:k
        for j in 1:k
            B[i,j] ~ Normal(0, 1)
        end
    end
    x = TV(undef, T, k)
    for t in 2:T
        # Debugging statements to check the types.
        x[t, :] ~ MvNormal(B * x[t-1, :], Σ_x)
        y[t, :] ~ MvNormal(x[t, :], Σ_y)
    end
end

m = SSM_2D(Y)
chain = sample(m, NUTS(0.8), 1000);
23

Can you say more about how it deviated from your expectations and what those expectations were?

24

Here is my MCMC implement result. means of parameters is not true parameter and
std is not estimated

parameters mean std mcse ess_bulk ess_tail rhat
F_x.L[1,1] 1.0000 0.0000 NaN NaN NaN NaN ⋯
F_x.L[2,1] 0.8824 0.0000 0.0000 NaN NaN NaN ⋯
F_x.L[2,2] 0.4706 0.0000 0.0000 NaN NaN NaN ⋯
σ_x[1] 0.1356 0.0000 0.0000 NaN NaN NaN ⋯
σ_x[2] 3.7133 0.0000 0.0000 NaN NaN NaN ⋯
F_y.L[1,1] 1.0000 0.0000 NaN NaN NaN NaN ⋯
F_y.L[2,1] 0.9602 0.0000 0.0000 NaN NaN NaN ⋯
F_y.L[2,2] 0.2793 0.0000 0.0000 NaN NaN NaN ⋯
σ_y[1] 5.6576 0.0000 0.0000 NaN NaN NaN ⋯
σ_y[2] 5.5809 0.0000 0.0000 NaN NaN NaN ⋯
B[1,1] 1.4806 0.0000 0.0000 NaN NaN NaN ⋯
B[1,2] -0.3448 0.0000 0.0000 NaN NaN NaN ⋯
B[2,1] -0.7867 0.0000 0.0000 NaN NaN NaN ⋯
B[2,2] 1.8545 0.0000 0.0000 NaN NaN NaN ⋯

25

I know it has been years but just revisiting this today and was able to use your suggestion here and everything works well with this toy example. Here it is in case anyone is interested:

using LinearAlgebra, PDMats, Random, Turing

Random.seed!(121)

# generate data  
σ = [1,2,3]  
Ω = [1 0.3 0.2;  
         0.3 1 0.1;  
         0.2 0.1 1]

Σ = Diagonal(σ) * Ω * Diagonal(σ)  
N = 100  
J = 3
y = rand(MvNormal(Σ), N)'

@model function correlation_chol(J, N, y)
	σ ~ filldist(truncated(Cauchy(0., 5.); lower = 0), J)
	F ~ LKJCholesky(J, 1.0)
    Σ_L = Diagonal(σ) * F.L
    Σ = PDMat(Cholesky(Σ_L + eps() * I))
    for i in 1:N
      y[i,:] ~ MvNormal(Σ)
    end
    # Compute Omega
    Ω = F.L * F.L'
    return (Σ=Σ, Ω=Ω)
end

model = correlation_chol(J, N, y)
chain = sample(model, NUTS(), 1000)

# check Omega results
gen_quants = generated_quantities(model, chain)
omega_matrices = [gq.Ω for gq in gen_quants]
mean(omega_matrices)

Resulting in:

image
image715×586 15.4 KB

and the following for Omega:

image

What’s inferred isn’t perfect but I think on-par with the Stan example I originally referenced. Of course increasing y length increases the accuracy in the results as well.

Thanks everyone for all the work you do for the open source community. Excited to start playing with Turing again.