Turing.jl: prophet NUTS sampling error no AdvancedHMC.DualValue(::Float64, ::Float64)

[bhillmann]

I am trying to reproduce the prophet model using Turing, but I am receiving an error when sampling with NUTS. I am fairly new to Julia so it could definitely be a user error. Currently, NUTS seems to fail to sample the trend component. I am able to sample using Metropilis Hastings. The current version of my code is below, but you can also find my environment on my repository https://github.com/bhillmann/juliacon2019/blob/master/turing_prophet_model.jl.

using Turing, StatsPlots, Random, CSV, DataFrames

Random.seed!(12)

# Reproducing the pymc3 implementation found here:
# https://www.ritchievink.com/blog/2018/10/09/build-facebooks-prophet-in-pymc3-bayesian-time-series-analyis-with-generalized-additive-models/

# downloaded from https://raw.githubusercontent.com/facebook/prophet/master/examples/example_wp_log_peyton_manning.csv
df = CSV.read("data/example_wp_log_peyton_manning.csv")

df[!, :yscale] = df[!, :y] ./ maximum(df[!, :y])

df[!, :t] = (df.ds - minimum(df.ds)) / (maximum(df.ds) - minimum(df.ds))

plot(df.ds, df.yscale)

nchangepoints = 10

# Proportion of history in which trend changepoints will be estimated.
changepoints_range = .8
# The standard deviation of the prior on the growth.
growth_prior_scale = 5
# The scale of the Laplace prior on the delta vector.
changepoints_prior_scale = .05

# Declare our Turing model.
@model trendmodel(t, y) = begin

    # num observations
    nobs = length(t)
    s =  collect(range(1, stop=nobs*changepoints_range, length=nchangepoints))
    A = [(t[i] >= s[j])*1 for i=1:nobs, j=1:nchangepoints]

    k ~ Normal(0.0, growth_prior_scale)

    # rate of change
    delta = Array{Real}(undef, nchangepoints)
    delta ~ [Laplace(0.0, changepoints_prior_scale)]

    # offset
    m ~ Normal(0.0, 5.0)

    gamma = -s.*delta

    g = (k .+ A*delta).*t + (m .+ A * gamma)

    # The number of observations.
    sd ~ Truncated(Cauchy(0.0, .5), 0.000001, Inf)
    for n in 1:nobs
        y[n] ~ Normal(g[n], sd)
    end
end

model = trendmodel(df.t, df.yscale)

chain = sample(model, NUTS(1500, 200, 0.65));

I receive the following error:

MethodError: no method matching AdvancedHMC.DualValue(::Float64, ::Float64)
Closest candidates are:
  AdvancedHMC.DualValue(::Tv<:AbstractFloat, !Matched::Tg<:AbstractArray{Tv<:AbstractFloat,1}) where {Tv<:AbstractFloat, Tg<:AbstractArray{Tv<:AbstractFloat,1}} at /home/benjamin/.julia/packages/AdvancedHMC/Ihu4Y/src/hamiltonian.jl:13
in top-level scope at base/none
in sample at Turing/rBXd7/src/inference/hmc.jl:228
in #sample#19 at Turing/rBXd7/src/inference/hmc.jl:279
in  at base/none
in #steps!#26 at Turing/rBXd7/src/inference/hmc.jl:421
in  at base/none
in #sample#11 at AdvancedHMC/Ihu4Y/src/sampler.jl:56
in phasepoint at AdvancedHMC/Ihu4Y/src/hamiltonian.jl:50
in ∂H∂θ at AdvancedHMC/Ihu4Y/src/hamiltonian.jl:20
Plots


REPL

[xukai92]

Also I remeber @torfjelde have some experience on implementing the model too. Can you take a look on this?

[xukai92]

I'm not sure if this is related to an known interface bug recently introduced. The bug is fixed not and we are wating for the version to be registered.

[xukai92]

(the order of comments above are supposed to be reversed)

[bhillmann]

So, should I update to the newest master and try it again?

[cpfiffer]

Might as well give it a shot, I think. You'll probably also need to be on the master for AdvancedHMC, though @xukai92 could speak to that a bit better.

[torfjelde]

Sorry, missed the mention; I'll dig up my code tomorrow :+1:

Glad to see you gave it a go @bhillmann!:)

[torfjelde]

My own (and old) implementation ran into the same issue, and it indeed seem to be the problem @xukai92 mentioned. Using Turing#master it works now:)

I ran the below code

@model trendmodel(t, y, ::Type{TV}=Vector{Float64}) where TV = begin  # <= CHANGED

    # num observations
    nobs = length(t)
    s =  collect(range(1, stop=nobs*changepoints_range, length=nchangepoints))
    A = [(t[i] >= s[j])*1 for i=1:nobs, j=1:nchangepoints]

    k ~ Normal(0.0, growth_prior_scale)

    # rate of change
    delta = TV(undef, nchangepoints)  # <= CHANGED
    delta ~ [Laplace(0.0, changepoints_prior_scale)]

    # offset
    m ~ Normal(0.0, 5.0)

    gamma = -s.*delta

    g = (k .+ A*delta).*t + (m .+ A * gamma)

    # The number of observations.
    sd ~ Truncated(Cauchy(0.0, .5), 0.000001, Inf)
    for n in 1:nobs
        y[n] ~ Normal(g[n], sd)
    end
end

Notice the # <= CHANGED comments. This functionality was added to Turing recently and makes the inference type-stable. Doing this you'll probably see something close to 10X speedup:) Also, on Turing#master, if you want to see the progress of the sampling add a progress = true to the sample call, i.e.

chain = sample(model, NUTS(1500, 200, 0.65), progress=true);

[xukai92]

Thanks for the help @torfjelde !

[garborg]

I tried to walk through the modified example in this issue, on Turing#master w/ AdvancedHMC#master on Julia 1.2.0, and the trendmodel call triggers a stackoverflow.

Here's a two-liner repro (w/ deps) if that helps lock down the issue: garborg/PNUTS.jl

[xukai92]

This is a know issue on Julia 1.2.0 (https://github.com/TuringLang/Turing.jl/issues/892).

The model should work with Julia 1.1.x.

[yebai]

Fixed by https://github.com/TuringLang/Turing.jl/pull/908https://github.com/TuringLang/Turing.jl/pull/908