# 2D state space model using Turing.jl
**URL:** https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486
**Category:** General Usage
**Created:** [October 2, 2023, 2:00am UTC](https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486 "2023-10-02T02:00:35Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@KotaAihara](https://discourse.julialang.org/u/KotaAihara)
#### Post date: [October 2, 2023, 2:00am UTC](https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486/1 "2023-10-02T02:00:35Z")
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.
```julia
#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?
---
### Author:  [@xukai92](https://discourse.julialang.org/u/xukai92)
#### Post date: [October 6, 2023, 2:16am UTC](https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486/2 "2023-10-06T02:16:07Z")
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](https://turinglang.org/v0.29/docs/using-turing/performancetips#ensure-that-types-in-your-model-can-be-inferred))
- 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:
```julia
@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
---
### Author:  [@KotaAihara](https://discourse.julialang.org/u/KotaAihara)
#### Post date: [October 6, 2023, 2:40am UTC](https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486/3 "2023-10-06T02:40:45Z")
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?
---
### Author:  [@xukai92](https://discourse.julialang.org/u/xukai92)
#### Post date: [October 6, 2023, 2:49am UTC](https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486/4 "2023-10-06T02:49:05Z")
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?
---
### Author:  [@KotaAihara](https://discourse.julialang.org/u/KotaAihara)
#### Post date: [October 11, 2023, 5:34am UTC](https://discourse.julialang.org/t/2d-state-space-model-using-turing-jl/104486/5 "2023-10-11T05:34:37Z")
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?
```julia
@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
```