Multivariate State Space model using Turing #2094

KotaAihara · 2023-10-03T11:32:53Z

KotaAihara
Oct 3, 2023

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?

torfjelde · 2023-10-03T13:18:12Z

torfjelde
Oct 3, 2023
Maintainer

x = tzeros(T, k) should not be used here; in fact, I don't think this is even defined anymore? Can you say which version of Julia and Turing.jl you're working with?

But to make this work, you can try:

@model function SSM_2D(y, ::Type{TV}=Matrix{Float64}) where {TV}
    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 = TV(undef, 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

9 replies

KotaAihara Oct 3, 2023
Author

I replaced it , but I got the same error. so I changed like this.

@model function SSM_2D(y, ::Type{TV}=Matrix{Float64}) where {TV}
    #y: observation data , n x m matrix
    T = size(y, 1) # T: number of observations
    k = size(y, 2) # k: number of variables
    #prior
    Σ_x ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    Σ_y ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    x = TV(undef, T, k)
    B = TV(undef, k, k)
    for i in 1:k
        for j in 1:k
            B[i,j] ~ Normal(0, 1)
        end
    end
    @show typeof(B)
    for t in 2:T
        μ_x = B * x[t-1, :]
        η_t ~ MvNormal(zeros(k), Σ_x)
        x[t, :] = μ_x + η_t
        ε_t ~ MvNormal(zeros(k), Σ_y)
        y[t, :] = x[t, :] + ε_t
    end
end

Unfortunately method error in ForwardDiff occurred .

torfjelde Oct 6, 2023
Maintainer

Can you post the first, say, 20 lines of the stack trace? I'm surprise the above results in issues.

torfjelde Oct 6, 2023
Maintainer

Oooh I see the issue. You're doing = for y.

Turing.jl only cares about lines involving ~, and so equalities are no handled. This means that you can't "infer" something on the right-hand side of =.

If you instead write:

        x[t, :] ~ MvNormal(B * x[t-1, :] + μ_x, Σ_x)
        y[t, :] ~ MvNormal(x[t, :] Σ_y)

in the loop, then it should work.

Also, in your example, Turing.jl will only see your η_t as a single variable, not a separate one for each timestep! This will completely break inference (I'm currently working on some features that will check for these issues before doing anything with the model, so these things should produce errors in the not-too-distant future).

torfjelde Oct 6, 2023
Maintainer

Ah, okay so there is in fact a bug in here.

I ran the following

using Turing, Random, Distributions, LinearAlgebra

T = 10

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
    Σ_x ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    Σ_y ~ InverseWishart(2, [1.0 0.0; 0.0 1.0])
    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.
        @info "" typeof(B * x[t-1, :])
        x[t, :] ~ MvNormal(B * x[t-1, :], Σ_x)
        @info "" typeof(x[t, :])
        y[t, :] ~ MvNormal(x[t, :], Σ_y)
    end
end

model = SSM_2D(Y)
# Execute the model to ensure that it's running without sampling.
model()

which resulted in

┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Any} (alias for Array{Any, 1})
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

Stacktrace:
  [1] macro expansion
    @ ./REPL[53]:18 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/DynamicPPL/m0PXI/src/compiler.jl:555 [inlined]
  [3] SSM_2D(__model__::DynamicPPL.Model{typeof(SSM_2D), (:y, Symbol("##arg#384")), (), (), Tuple{Matrix{Float64}, DataType}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, __context__::DynamicPPL.SamplingContext{DynamicPPL.SampleFromPrior, DynamicPPL.DefaultContext, TaskLocalRNG}, y::Matrix{Float64}, arg#384::Type{Matrix{Float64}})
    @ Main ./REPL[53]:14
  [4] _evaluate!!(model::DynamicPPL.Model{typeof(SSM_2D), (:y, Symbol("##arg#384")), (), (), Tuple{Matrix{Float64}, DataType}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, context::DynamicPPL.SamplingContext{DynamicPPL.SampleFromPrior, DynamicPPL.DefaultContext, TaskLocalRNG})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/m0PXI/src/model.jl:963
  [5] evaluate_threadunsafe!!
    @ ~/.julia/packages/DynamicPPL/m0PXI/src/model.jl:936 [inlined]
  [6] evaluate!!
    @ ~/.julia/packages/DynamicPPL/m0PXI/src/model.jl:889 [inlined]
  [7] evaluate!! (repeats 2 times)
    @ ~/.julia/packages/DynamicPPL/m0PXI/src/model.jl:900 [inlined]
  [8] evaluate!!
    @ ~/.julia/packages/DynamicPPL/m0PXI/src/model.jl:908 [inlined]
  [9] (::DynamicPPL.Model{typeof(SSM_2D), (:y, Symbol("##arg#384")), (), (), Tuple{Matrix{Float64}, DataType}, Tuple{}, DynamicPPL.DefaultContext})()
    @ DynamicPPL ~/.julia/packages/DynamicPPL/m0PXI/src/model.jl:860
 [10] top-level scope
    @ REPL[55]:1

This is indeed a bug. You can fix this by adding BangBang.jl to your enviroment, adding the following snippet:

using BangBang

function BangBang.possible(
    ::typeof(BangBang._setindex!), ::C, ::T, ::Integer, ::Colon
    ) where {C<:AbstractMatrix,T<:AbstractVector}
    return BangBang.implements(setindex!, C) &&
        promote_type(eltype(C), eltype(T)) <: eltype(C)
end

Once this has been done, I now get:

julia> model()
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(B * x[t - 1, :]) = Vector{Float64} (alias for Array{Float64, 1})
┌ Info: 
└   typeof(x[t, :]) = Vector{Float64} (alias for Array{Float64, 1})

Removing the @info statements:

julia> @model function SSM_2D(y, ::Type{TV}=Matrix{Float64}) where {TV}
           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 = 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

SSM_2D (generic function with 4 methods)

julia> model = SSM_2D(Y);

julia> chain = sample(model, NUTS(0.65), 1000);
┌ Info: Found initial step size
└   ϵ = 0.00625
Sampling 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:18

and it works!

I'll make a PR to fix this though; thanks for bringing this to our attention:)

KotaAihara Oct 10, 2023
Author

I appreciate your assistance. Occasionally, I have encountered a message regarding a failure in Cholesky Decomposition, so I adjusted the prior of the Variance-Covariance matrix. However, the posterior of the coefficient matrix differs from what I inputted. Do you have any ideas?

#seed
Random.seed!(1234)

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.65), 1000);

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Multivariate State Space model using Turing #2094

{{title}}

Replies: 1 comment 9 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Multivariate State Space model using Turing #2094

KotaAihara Oct 3, 2023

Replies: 1 comment · 9 replies

torfjelde Oct 3, 2023 Maintainer

KotaAihara Oct 3, 2023 Author

torfjelde Oct 6, 2023 Maintainer

torfjelde Oct 6, 2023 Maintainer

torfjelde Oct 6, 2023 Maintainer

KotaAihara Oct 10, 2023 Author

KotaAihara
Oct 3, 2023

Replies: 1 comment 9 replies

torfjelde
Oct 3, 2023
Maintainer

KotaAihara Oct 3, 2023
Author

torfjelde Oct 6, 2023
Maintainer

torfjelde Oct 6, 2023
Maintainer

torfjelde Oct 6, 2023
Maintainer

KotaAihara Oct 10, 2023
Author