# Issues with HMC/NUTS samplers

**URL:** <https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900>\
**Category:** Probabilistic Programming\
**Tags:** turing, dynamichmc\
**Created:** [September 29, 2021, 5:15am UTC](https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900 "2021-09-29T05:15:20Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![tchebycheff](https://avatars.discourse-cdn.com/v4/letter/t/779978/32.png) [@tchebycheff](https://discourse.julialang.org/u/tchebycheff)\
**Post date:** [September 29, 2021, 5:15am UTC](https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900/1 "2021-09-29T05:15:20Z")

</div>

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

```julia
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.

```julia
@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

```

---

<div class="post-metadata">

**Author:** ![sami1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sami1/32/23555_2.png) [@sami1](https://discourse.julialang.org/u/sami1)\
**Post date:** [September 29, 2021, 10:42am UTC](https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900/2 "2021-09-29T10:42:35Z")

</div>

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 :

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

```

and pass your parameters through it.

---

<div class="post-metadata">

**Author:** ![tchebycheff](https://avatars.discourse-cdn.com/v4/letter/t/779978/32.png) [@tchebycheff](https://discourse.julialang.org/u/tchebycheff)\
**Post date:** [September 29, 2021, 1:50pm UTC](https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900/3 "2021-09-29T13:50:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [September 29, 2021, 2:01pm UTC](https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900/4 "2021-09-29T14:01:41Z")

</div>

> [@tchebycheff](#):
>
> It works now.

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

---

<div class="post-metadata">

**Author:** ![tchebycheff](https://avatars.discourse-cdn.com/v4/letter/t/779978/32.png) [@tchebycheff](https://discourse.julialang.org/u/tchebycheff)\
**Post date:** [September 29, 2021, 3:53pm UTC](https://discourse.julialang.org/t/issues-with-hmc-nuts-samplers/68900/5 "2021-09-29T15:53:41Z")

</div>

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

```julia
@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

```