Skip to content

Instantly share code, notes, and snippets.

@JasonPekos
Created May 26, 2024 17:06
Show Gist options
  • Save JasonPekos/182aa1c880676faca92c28c744ddf084 to your computer and use it in GitHub Desktop.
Save JasonPekos/182aa1c880676faca92c28c744ddf084 to your computer and use it in GitHub Desktop.
using HiddenMarkovModels, Turing, LinearAlgebra
using FillArrays, LogExpFunctions
# This is a toy example taken from the Turing tutorial on HMMs:
# https://turinglang.org/docs/tutorials/04-hidden-markov-model/
# Three possible states, emission are just normally distributed around the latent mean.
obs = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0,
3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,];
K = 3; # Number of states.
@model function hmm_working(y, K)
m ~ MvNormal([1, 2, 3], I)
T ~ filldist(Dirichlet(Fill(1/K, K)), K)
hmm = HMM(softmax(ones(K)), copy(T'), [Normal(m[i], 0.1) for i in 1:K])
Turing.@addlogprob! logdensityof(hmm, y)
end
@model function hmm_not_working(y, K)
m ~ MvNormal([1, 2, 3], I)
T ~ filldist(Dirichlet(Fill(1/K, K)), K)
# the difference is the T' vs copy(T') on this line:
hmm = HMM(softmax(ones(K)), T', [Normal(m[i], 0.1) for i in 1:K])
Turing.@addlogprob! logdensityof(hmm, y)
end
sample(hmm_working(obs, K), NUTS(), 100) # this sample call works
sample(hmm_not_working(obs, K), NUTS(), 100) # this one doesn't
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
HiddenMarkovModels = "84ca31d5-effc-45e0-bfda-5a68cd981f47"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
ERROR: MethodError: Cannot `convert` an object of type Matrix{Float64} to an object of type Adjoint{Float64, Matrix{Float64}}
Closest candidates are:
convert(::Type{Adjoint{T, S}}, ::Adjoint) where {T, S}
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/adjtrans.jl:337
convert(::Type{T}, ::T) where T
@ Base Base.jl:84
convert(::Type{T}, ::Factorization) where T<:AbstractArray
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:108
...
Stacktrace:
[1] HMM(init::Vector{Float64}, trans::Adjoint{Float64, Matrix{Float64}}, dists::Vector{Normal{Float64}})
@ HiddenMarkovModels ~/.julia/packages/HiddenMarkovModels/bhz3w/src/types/hmm.jl:23
[2] hmm_not_working(__model__::DynamicPPL.Model{…}, __varinfo__::DynamicPPL.ThreadSafeVarInfo{…}, __context__::DynamicPPL.SamplingContext{…}, y::Vector{…}, K::Int64)
@ Main ~/Documents/julia-scratch-dev/hmms-turing/min-rep.jl:28
[3] _evaluate!!(model::DynamicPPL.Model{…}, varinfo::DynamicPPL.ThreadSafeVarInfo{…}, context::DynamicPPL.SamplingContext{…})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:963
[4] evaluate_threadsafe!!(model::DynamicPPL.Model{…}, varinfo::DynamicPPL.UntypedVarInfo{…}, context::DynamicPPL.SamplingContext{…})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:952
[5] evaluate!!(model::DynamicPPL.Model{…}, varinfo::DynamicPPL.UntypedVarInfo{…}, context::DynamicPPL.SamplingContext{…})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:887
[6] evaluate!!(model::DynamicPPL.Model{…}, rng::Random.TaskLocalRNG, varinfo::DynamicPPL.UntypedVarInfo{…}, sampler::DynamicPPL.SampleFromUniform, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:900
[7] (::DynamicPPL.Model{…})(::Random.TaskLocalRNG, ::Vararg{…})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:860
[8] DynamicPPL.VarInfo(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.SampleFromUniform, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/varinfo.jl:129
[9] default_varinfo(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…}, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/sampler.jl:80
[10] default_varinfo(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/sampler.jl:71
[11] step(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, spl::DynamicPPL.Sampler{…}; initial_params::Nothing, kwargs::@Kwargs{…})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/sampler.jl:103
[12] macro expansion
@ ~/.julia/packages/AbstractMCMC/YrmkI/src/sample.jl:130 [inlined]
[13] macro expansion
@ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
[14] macro expansion
@ ~/.julia/packages/AbstractMCMC/YrmkI/src/logging.jl:9 [inlined]
[15] mcmcsample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, initial_state::Nothing, kwargs::@Kwargs{…})
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/YrmkI/src/sample.jl:120
[16] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…}, N::Int64; chain_type::Type, resume_from::Nothing, initial_state::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::@Kwargs{})
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/hmc.jl:117
[17] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/hmc.jl:86
[18] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, alg::NUTS{…}, N::Int64; kwargs::@Kwargs{})
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:263
[19] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{…}, alg::NUTS{…}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:256
[20] sample(model::DynamicPPL.Model{…}, alg::NUTS{…}, N::Int64; kwargs::@Kwargs{})
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:253
[21] sample(model::DynamicPPL.Model{…}, alg::NUTS{…}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:247
[22] top-level scope
@ ~/Documents/julia-scratch-dev/hmms-turing/min-rep.jl:34
Some type information was truncated. Use `show(err)` to see complete types.
1-element ExceptionStack:
LoadError: MethodError: Cannot `convert` an object of type Matrix{Float64} to an object of type Adjoint{Float64, Matrix{Float64}}
Closest candidates are:
convert(::Type{Adjoint{T, S}}, ::Adjoint) where {T, S}
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/adjtrans.jl:337
convert(::Type{T}, ::T) where T
@ Base Base.jl:84
convert(::Type{T}, ::Factorization) where T<:AbstractArray
@ LinearAlgebra ~/.julia/juliaup/julia-1.10.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:108
...
Stacktrace:
[1] HMM(init::Vector{Float64}, trans::Adjoint{Float64, Matrix{Float64}}, dists::Vector{Normal{Float64}})
@ HiddenMarkovModels ~/.julia/packages/HiddenMarkovModels/bhz3w/src/types/hmm.jl:23
[2] hmm_not_working(__model__::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random.TaskLocalRNG}, y::Vector{Float64}, K::Int64)
@ Main ~/Documents/julia-scratch-dev/hmms-turing/min-rep.jl:28
[3] _evaluate!!(model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}}, context::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random.TaskLocalRNG})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:963
[4] evaluate_threadsafe!!(model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, 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.SampleFromUniform, DynamicPPL.DefaultContext, Random.TaskLocalRNG})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:952
[5] evaluate!!(model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, 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.SampleFromUniform, DynamicPPL.DefaultContext, Random.TaskLocalRNG})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:887
[6] evaluate!!(model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, rng::Random.TaskLocalRNG, varinfo::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, sampler::DynamicPPL.SampleFromUniform, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:900
[7] (::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext})(::Random.TaskLocalRNG, ::Vararg{Any})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/model.jl:860
[8] DynamicPPL.VarInfo(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.SampleFromUniform, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/varinfo.jl:129
[9] default_varinfo(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, context::DynamicPPL.DefaultContext)
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/sampler.jl:80
[10] default_varinfo(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/sampler.jl:71
[11] step(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}; initial_params::Nothing, kwargs::@Kwargs{nadapts::Int64})
@ DynamicPPL ~/.julia/packages/DynamicPPL/E4kDs/src/sampler.jl:103
[12] macro expansion
@ ~/.julia/packages/AbstractMCMC/YrmkI/src/sample.jl:130 [inlined]
[13] macro expansion
@ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
[14] macro expansion
@ ~/.julia/packages/AbstractMCMC/YrmkI/src/logging.jl:9 [inlined]
[15] mcmcsample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, initial_state::Nothing, kwargs::@Kwargs{nadapts::Int64})
@ AbstractMCMC ~/.julia/packages/AbstractMCMC/YrmkI/src/sample.jl:120
[16] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; chain_type::Type, resume_from::Nothing, initial_state::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::@Kwargs{})
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/hmc.jl:117
[17] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/hmc.jl:86
[18] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64; kwargs::@Kwargs{})
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:263
[19] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:256
[20] sample(model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64; kwargs::@Kwargs{})
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:253
[21] sample(model::DynamicPPL.Model{typeof(hmm_not_working), (:y, :K), (), (), Tuple{Vector{Float64}, Int64}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{AutoForwardDiff{nothing, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64)
@ Turing.Inference ~/.julia/packages/Turing/cH3wV/src/mcmc/Inference.jl:247
[22] top-level scope
@ ~/Documents/julia-scratch-dev/hmms-turing/min-rep.jl:34
[23] eval
@ ./boot.jl:385 [inlined]
[24] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:2076
[25] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
@ Base ./essentials.jl:892
[26] invokelatest(::Any, ::Any, ::Vararg{Any})
@ Base ./essentials.jl:889
[27] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:271
[28] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:181
[29] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/repl.jl:276
[30] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:179
[31] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/repl.jl:38
[32] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:150
[33] with_logstate(f::Function, logstate::Any)
@ Base.CoreLogging ./logging.jl:515
[34] with_logger
@ ./logging.jl:627 [inlined]
[35] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:263
[36] #invokelatest#2
@ ./essentials.jl:892 [inlined]
[37] invokelatest(::Any)
@ Base ./essentials.jl:889
[38] (::VSCodeServer.var"#64#65")()
@ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.79.2/scripts/packages/VSCodeServer/src/eval.jl:34
in expression starting at /Users/pekos/Documents/julia-scratch-dev/hmms-turing/min-rep.jl:34
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment