Issues with HMC/NUTS samplers

Specific Domains Probabilistic programming
,
1

Hi - I am trying to run MCMC on this model with Turing.
Only SMC sampler works. With HMC, NUTS, PG, DynamicHMC I get the error

TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 10}
    Stacktrace:

I experimented with different parameters for samplers. No cigar. Any suggestions/advice is much appreciated. Thanks.

Model:
y_t = s_t * n_t
n_t is iid N(0,1) and s_t = s0 * sqrt(M1_t * M2_t …Mk_t)
Mj_t’s are independent (2-state) Markov chains, taking non negative values.

@model function func(y, kbar)
    N = length(y)
    S = Vector{Vector{Int64}}(undef, N)
    A = Vector{Matrix{Float64}}(undef, kbar)

    b ~ Uniform(1.0, 50.)
    m0 ~ Uniform(1.0, 1.9999)
    γk ~ Uniform(0.001, 0.999)
    σ0 ~ Uniform(0.0001, 5.)

    M = [m0, 2. - m0]

    for i = 1:kbar
        γ = (1. - (1. - γk)^(b^(i - kbar))) / 2.
        A[i] = [1.0 - γ γ; γ 1.0 - γ]
    end

    S[1] ~ Product(Categorical.(ones(Int64, kbar) * 2))
    σ = σ0 * sqrt(prod([M[u] for u in S[1]]))
    y[1] ~ Normal(0.0, σ)

    for j = 2:N
        S[j] ~ Product(Categorical.([A[i][S[j - 1][i],:] for i = 1:kbar]))
        σ = σ0 * sqrt(prod([M[u] for u in S[j - 1]]))
        y[j] ~ Normal(0, σ)
    end
    return y, S
end
2

Hi,
This is a common type issue with autodifferentiation (that’s why you didn’t encounter it with ohter samplers).
Your vector A is pre-allocated with a Matrix{Float64} as element type, but later in the for loop you update it with [1.0 - γ γ; γ 1.0 - γ] which is of type Matrix{ForwardDiff.Dual}.
A quick fix is to avoid specifying the element type when you initialize arrays.

If you need type specification (e.g. for multiple dispatch), you can use a function like :

check_type(x) = x
check_type(x::T) where T <:ForwardDiff.Dual = x.value

and pass your parameters through it.

2 Likes
3

Thanks Sami. That helped.
There was also an issue with using the Categorical rv S as index to a vector (M in the code).
It works now.

4

For the benefit of the rest of us, can you show what the working code looks like?

5

Hi -
This is the corrected version. There is probably plenty of room for optimisation. Would appreciate any suggestions.

@model function msm(y, kbar)
    N = length(y)
    A = Vector{Matrix}(undef, kbar)
    S = Vector(undef, kbar)

    b ~ Uniform(1.0, 50.)
    m0 ~ Uniform(1.0, 1.9999)
    γ1 ~ Uniform(0.001, 0.999)
    σ0 ~ Uniform(0.0001, 5.)

    M = [m0, 2.0 - m0]

    for i = 1:kbar
        γ = (1. - (1. - γ1)^(b^(i - 1))) / 2.
        A[i] = [1.0 - γ γ; γ 1.0 - γ]
    end
    S ~ product_distribution([Categorical([0.5, 0.5]) for i = 1:kbar])
    σ = σ0 * sqrt(prod([M[u] for u in S]))
    y[1] ~ Normal(0.0, σ)

   for j = 2:N
        S ~ product_distribution([Categorical(A[i][S[i],:]) for i = 1:kbar])
        σ = σ0 * sqrt(prod([M[u] for u in S]))
        y[j] ~ Normal(0.0, σ)
    end
end