# 2D state space model using Turing.jl

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

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})
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?

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

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?

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?

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

``````