使用 Turing.jl 对 AR(2) 过程进行采样时出现类型错误

编程问答 2022-05-09

问题描述

我正在尝试创建一个简单的 AR(2) 流程 Y 其中：

Y[t] = ϵ[t] + b1 * y[t - 1] + b2 * y[t - 2]
ϵ[t] ~ normal(0,σ)

、b1 和 b2 是从一些先验分布中提取的参数。

代码如下：

using Statistics
using Turing,distributions

# This is supposed to be AR(2) model
# Y[t] = ε[t] + b1 * Y[t-1] + b2 * Y[t-2]
@model function AR2(y,b1_mean,b2_mean,::Type{T}=Float64; σ=.3,σ_b=.01) where T <: Real
    N = length(y)
    
    b1 ~ normal(b1_mean,σ_b)
    b2 ~ normal(b2_mean,σ_b)
    ϵ = Vector{T}(undef,N)
    
    y[1] ~ normal(0,σ)
    y[2] ~ normal(0,σ)
    
    for k ∈ 3:N
        ϵ[k] ~ normal(0,σ)
        # Maybe I shouldn't be assigning here?
        y[k] = ϵ[k] + b1 * y[k - 1] + b2 * y[k - 2] # this is LINE 19
    end
    
    y
end

# These are the parameters
b1,b2 = .3,.6

println("Constructing random process with KNowN parameters: b1=$b1,b2=$b2")
AR2_process = AR2(fill(missing,500),b1,b2)()
AR2_process_orig = copy(AR2_process)

@show typeof(AR2_process)

println("Estimating parameters...")
model_ar = AR2(AR2_process,b2)
chain_ar = sample(model_ar,HMC(0.001,10),200)

@show mean(chain_ar[:b1]),mean(chain_ar[:b2])
@show all(AR2_process .≈ AR2_process_orig)

println("Now try to estimate Float64")
model_ar = AR2(Float64.(AR2_process),mean(chain_ar[:b2])

输出是这样的：

Constructing random process with KNowN parameters: b1=0.3,b2=0.6
typeof(AR2_process) = Vector{Union{Missing,Float64}}
Estimating parameters...
Sampling 100%|█████████████████| Time: 0:00:47
(mean(chain_ar[:b1]),mean(chain_ar[:b2])) = (0.29404359023752596,0.5912394358365198)
all(AR2_process .≈ AR2_process_orig) = true
Now try to estimate Float64
Sampling 100%|█████████████████| Time: 0:00:02
ERROR: LoadError: TypeError: in typeassert,expected Float64,got a value of type ForwardDiff.Dual{nothing,Float64,10}
Stacktrace:
  [1] setindex!(A::Vector{Float64},x::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,:b2,:ϵ),Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:b1,Tuple{}},Int64},Vector{normal{Float64}},Vector{AbstractPPL.VarName{:b1,Tuple{}}},Vector{Float64},Vector{Set{DynamicPPL.Selector}}},DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:b2,Vector{AbstractPPL.VarName{:b2,DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ϵ,Tuple{Tuple{Int64}}},Vector{AbstractPPL.VarName{:ϵ,Tuple{Tuple{Int64}}}},Vector{Set{DynamicPPL.Selector}}}}},Float64},DynamicPPL.Model{var"#2#3",(:y,:b1_mean,:b2_mean,:T,:σ,:σ_b),(:T,(),Tuple{Vector{Float64},Type{Float64},Tuple{Type{Float64},Float64}},DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},AdvancedHMC.UnitEuclideanMetric}},DynamicPPL.DefaultContext},10},i1::Int64)
    @ Base ./array.jl:839
  [2] #2
    @ ~/test/turing_question.jl:19 [inlined]
  [3] (::var"#2#3")(__rng__::Random._GLOBAL_RNG,__model__::DynamicPPL.Model{var"#2#3",__varinfo__::DynamicPPL.TypedVarInfo{NamedTuple{(:b1,Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,10}},ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,__sampler__::DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},__context__::DynamicPPL.DefaultContext,y::Vector{Float64},b1_mean::Float64,b2_mean::Float64,#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,σ::Float64,σ_b::Float64)
    @ Main ./none:0
  [4] macro expansion
    @ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:0 [inlined]
  [5] _evaluate
    @ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:154 [inlined]
  [6] evaluate_threadunsafe
    @ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:127 [inlined]
  [7] (::DynamicPPL.Model{var"#2#3",Float64}})(rng::Random._GLOBAL_RNG,varinfo::DynamicPPL.TypedVarInfo{NamedTuple{(:b1,sampler::DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},context::DynamicPPL.DefaultContext)
    @ DynamicPPL ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:92
  [8] Model
    @ ~/.julia/packages/DynamicPPL/SgzCy/src/model.jl:98 [inlined]
  [9] f
    @ ~/.julia/packages/Turing/28kgo/src/core/ad.jl:111 [inlined]
 [10] chunk_mode_gradient!(result::Vector{Float64},f::Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,x::Vector{Float64},cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,10,10}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:150
 [11] gradient!
    @ ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:39 [inlined]
 [12] gradient!(result::Vector{Float64},10}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:35
 [13] gradient_logp(::Turing.Core.ForwardDiffAD{40},θ::Vector{Float64},vi::DynamicPPL.TypedVarInfo{NamedTuple{(:b1,model::DynamicPPL.Model{var"#2#3",ctx::DynamicPPL.DefaultContext)
    @ Turing.Core ~/.julia/packages/Turing/28kgo/src/core/ad.jl:121
 [14] gradient_logp (repeats 2 times)
    @ ~/.julia/packages/Turing/28kgo/src/core/ad.jl:83 [inlined]
 [15] ∂logπ∂θ
    @ ~/.julia/packages/Turing/28kgo/src/inference/hmc.jl:433 [inlined]
 [16] ∂H∂θ
    @ ~/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:31 [inlined]
 [17] phasepoint(h::AdvancedHMC.Hamiltonian{AdvancedHMC.UnitEuclideanMetric{Float64,Tuple{Int64}},Turing.Inference.var"#logπ#52"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,Float64}}},Turing.Inference.var"#∂logπ∂θ#51"{DynamicPPL.TypedVarInfo{NamedTuple{(:b1,Float64}}}},r::Vector{Float64})
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:69
 [18] phasepoint
    @ ~/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:139 [inlined]
 [19] initialstep(rng::Random._GLOBAL_RNG,spl::DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Float64}; init_params::nothing,nadapts::Int64,kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}})
    @ Turing.Inference ~/.julia/packages/Turing/28kgo/src/inference/hmc.jl:167
 [20] initialstep(rng::Random._GLOBAL_RNG,Float64})
    @ Turing.Inference ~/.julia/packages/Turing/28kgo/src/inference/hmc.jl:153
 [21] step(rng::Random._GLOBAL_RNG,AdvancedHMC.UnitEuclideanMetric}}; resume_from::nothing,Tuple{}}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/SgzCy/src/sampler.jl:91
 [22] step
    @ ~/.julia/packages/DynamicPPL/SgzCy/src/sampler.jl:66 [inlined]
 [23] macro expansion
    @ ~/.julia/packages/AbstractMCMC/ByHEr/src/sample.jl:123 [inlined]
...

（因为帖子超过了 30000 个字符的限制，我不得不删掉错误信息）

所以它是说第 19 行：

y[k] = ϵ[k] + b1 * y[k - 1] + b2 * y[k - 2]

...试图将“ForwardDiff.Dual{nothing,10} 类型的值”插入到 y 中，但 y 是 Vector{Float64} 类型，因为我明确地将其强制转换为 {{ 1}} 与 Float64，因此类型不匹配导致错误。这有点道理。

问题

为什么上一个运行有效？为什么从 Float64.(AR2_process) 采样（其中 AR2(AR2_process,b2) 是 AR2_process）不会引发任何错误？它不会在这里也尝试将 Vector{Union{Missing,Float64}} 分配给 ForwardDiff.Dual 吗？

潜在的解决方案

也许我不应该将内容分配给y[k]，因为观测值应该从某个分布抽取，所以它可能应该是{{1 }} 因为我正在添加从正态分布中提取的值，但我不太确定该正态分布的参数应该是什么...

解决方法

我不知道为什么您的代码会以这种方式失败，但我可以解释并解决您的实际问题。

是的，您必须始终使用波浪线进行观察；但关键是您的观察陈述的右侧不是分布，而是确定性的。如果您实际上不需要采样/估计 epsilon，则可以改用以下模型：

@model function AR2(y,b1_mean,b2_mean; σ=.3,σ_b=.01) where T <: Real
    N = length(y)
           
    b1 ~ Normal(b1_mean,σ_b)
    b2 ~ Normal(b2_mean,σ_b)
           
    y[1] ~ Normal(0,σ)
    y[2] ~ Normal(0,σ)
           
    for k ∈ 3:N
        y[k] ~ Normal(b1 * y[k - 1] + b2 * y[k - 2],σ) # this is LINE 19
    end
           
    y
end

这只是因为我们知道 Normal 的属性。跟踪随机变量的确定性函数的问题仍然是an open issue。

或者，您尝试新的热门内容并编写基于 submodel 的实现。

juliastochastic-process