2D state space model using Turing.jl

General Usage
1

I’m making 2d state space model by using Turing.jl

Here is my code. x[t] is state equation and y[t] is observation equation.

#Multivariate SSM
#sample model 

# model 
using Turing 
using Random, Distributions
using Plots
using StaticArrays,LinearAlgebra
using Libtask

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

#plot data
plot(Y[:, 1], label="Y1", title="Multivariate State Space Model", xlabel="t", ylabel="Y", legend=:bottomright)
plot!(Y[:, 2], label="Y2")
plot!(X[:, 1], label="X1")
plot!(X[:, 2], label="X2")


@model function SSM_2D(y)
    T = size(y, 1)
    k = size(y, 2)
    #prior
    Σ_x ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    Σ_y ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    B = zeros(k, k  )
    for i in 1:k
        for j in 1:k
            B[i,j] ~ Normal(0, 1)
        end
    end
    x = tzeros(T, k)
    for t in 2:T
        μ_x = B * x[t-1, :]
        x[t, :] ~ MvNormal(μ_x, Σ_x)
        x_t = x[t, :]
        y[t, :] ~ MvNormal(x_t, Σ_y)
    end
end

#MCMC sampling
model = SSM_2D(Y)
chain = sample(model, NUTS(0.65), 1000, progress=true)

but I got this error.

MethodError: no method matching MvNormal(::Vector{Any}, ::Matrix{Float64})

Closest candidates are:
MvNormal(::Tracker.TrackedVector{<:Real}, ::Matrix{<:Real})
@ DistributionsADTrackerExt ~/.julia/packages/DistributionsAD/rvsL2/ext/DistributionsADTrackerExt.jl:441
MvNormal(::AbstractVector{<:Real}, ::AbstractMatrix{<:Real})
@ Distributions ~/.julia/packages/Distributions/fYgbJ/src/multivariate/mvnormal.jl:201

I suppose it happened in for loop. Could you teach me how to fix my error?

2

the direct cause is that x = tzeros(T, k) creates an array of type Any that MvNormal does not like and you are also trying to assign a vector to a matrix row via x[t, :] ~ MvNormal(μ_x, Σ_x) which triggers type conversion and results Any

the fix is as fllow:

  • updating the model signature to function SSM_2D(y, ::Type{FT}=Float64) where FT and using x = tzeros(FT, T, k) (ref: Performance Tips)
  • this will makes sure the array is correctly typed
  • also, you don’t need to use tzeros actually if you are not using particle-based samplers, which seem to be not necessary because you are only dealing with continuous variables and using HMC (NUTS) so x = zeros(FT, T, k) should be fine too
  • make matrices column-major (Julia’s preference)

complete model def below:

@model function SSM_2D(y, ::Type{FT}=Float64) where {FT}
    T = size(y, 1)
    k = size(y, 2)
    #prior
    Σ_x ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    Σ_y ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    B = zeros(FT, k, k)
    for i in 1:k
        for j in 1:k
            B[i,j] ~ Normal(0, 1)
        end
    end
    x = zeros(FT, k, T)
    for t in 2:T
        μ_x = B * x[:, t-1]
        x[:, t] ~ MvNormal(μ_x, Σ_x)
        x_t = x[:, t]
        y[t, :] ~ MvNormal(x_t, Σ_y)
    end
end
  • i didn’t change that for y but feel free to do the same column-major update, along with the data generation
  • i updated B = zeros(FT, k, k) as well for the type thing
3

Thank you ! it worked.
But I got this error

PosDefException: matrix is not positive definite; Cholesky factorization failed.

Do you have some idea to constrain parameters?

4

interesting.
i was able to run the code with T = 5 in the data simulation part (it was missing in the OP) but looks like if i use a larger T it fails.
would it be related to parameters in the ground truth that makes it hard or something with your prior?

5

I am grateful for your assistance. There have been instances when I received a message concerning a failure in Cholesky Decomposition, prompting me to modify the prior of the Variance-Covariance matrix. Nonetheless, the posterior of the coefficient matrix deviates from what I had inputted. Would you have any insights?

@model function SSM_2D(y, ::Type{FT}=Float64) where {FT}
    T = size(y, 1)
    k = size(y, 2)
    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 = zeros(FT, k, k)
    for i in 1:k
        for j in 1:k
            B[i,j] ~ Normal(0, 1)
        end
    end
    x = zeros(FT, k, T)
    for t in 2:T
        μ_x = B * x[:, t-1]
        x[:, t] ~ MvNormal(μ_x, Σ_x)
        x_t = x[:, t]
        y[t, :] ~ MvNormal(x_t, Σ_y)
    end
end