DynamicHMC: Reached maximum number of iterations

General Usage
1

This is likely a question for Tamas. @goedman and I have encountered an error the following error with DynamicHMC: ERROR: Reached maximum number of iterations while bisecting interval for ϵ. The error occurs intermittently, about 1 in 5 runs. Here is the full error message:

ERROR: Reached maximum number of iterations while bisecting interval for ϵ.
Stacktrace:
 [1] bisect_stepsize(::DynamicHMC.InitialStepsizeSearch, ::getfield(DynamicHMC, Symbol("##3#4")){DynamicHMC.Hamiltonian{LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6},1}}},GaussianKE{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.LowerTriangular{Float64,LinearAlgebra.Diagonal{Float64,Array{Float64,1}}}}},DynamicHMC.PhasePoint{Array{Float64,1},LogDensityProblems.ValueGradient{Float64,Array{Float64,1}}},Float64}, ::Float64, ::Float64, ::Float64, ::Float64) at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/stepsize.jl:122
 [2] find_initial_stepsize at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/stepsize.jl:146 [inlined]
 [3] find_initial_stepsize(::DynamicHMC.InitialStepsizeSearch, ::DynamicHMC.Hamiltonian{LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6},1}}},GaussianKE{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.LowerTriangular{Float64,LinearAlgebra.Diagonal{Float64,Array{Float64,1}}}}}, ::DynamicHMC.PhasePoint{Array{Float64,1},LogDensityProblems.ValueGradient{Float64,Array{Float64,1}}}) at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/stepsize.jl:151
 [4] #NUTS_init#18(::Array{Float64,1}, ::GaussianKE{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.LowerTriangular{Float64,LinearAlgebra.Diagonal{Float64,Array{Float64,1}}}}, ::Array{Float64,1}, ::Int64, ::DynamicHMC.InitialStepsizeSearch, ::ReportIO{Base.TTY}, ::Function, ::MersenneTwister, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6},1}}}) at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/sampler.jl:141
 [5] NUTS_init(::MersenneTwister, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6},1}}}) at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/sampler.jl:138
 [6] #NUTS_init_tune_mcmc#20(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::MersenneTwister, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6},1}}}, ::Int64) at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/sampler.jl:284
 [7] NUTS_init_tune_mcmc at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/sampler.jl:284 [inlined]
 [8] #NUTS_init_tune_mcmc#21 at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/sampler.jl:294 [inlined]
 [9] NUTS_init_tune_mcmc(::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:v, :A, :k, :tau),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64}}}},LBAProb{NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}}}},Float64},Float64,6},1}}}, ::Int64) at /home/dfish/.julia/packages/DynamicHMC/YX6JA/src/sampler.jl:294
 [10] sampleDHMC(::NamedTuple{(:choice, :rt),Tuple{Array{Int64,1},Array{Float64,1}}}, ::Int64, ::Int64, ::Int64) at /home/dfish/.julia/dev/MCMCBenchmarks/src/temp.jl:159
 [11] top-level scope at none:0

We were wondering if you have any advice for dealing with this situation. One solution that comes immediately to mind is to increase the number of iterations. However, that may not be the most efficient solution. Here is a MWE (sorry the model definitions are long):

using Distributions,Parameters,DynamicHMC,LogDensityProblems,TransformVariables
using Random
import Distributions: pdf,logpdf,rand
export LBA,pdf,logpdf,rand

############################################################################################
#                                    Model functions
############################################################################################

mutable struct LBA{T1,T2,T3,T4} <: ContinuousUnivariateDistribution
    ν::T1
    A::T2
    k::T3
    τ::T4
    σ::Float64
end

Base.broadcastable(x::LBA)=Ref(x)

LBA(;τ,A,k,ν,σ=1.0) = LBA(ν,A,k,τ,σ)

function selectWinner(dt)
    if any(x->x >0,dt)
        mi,mv = 0,Inf
        for (i,t) in enumerate(dt)
            if (t > 0) && (t < mv)
                mi = i
                mv = t
            end
        end
    else
        return 1,-1.0
    end
    return mi,mv
end

function sampleDriftRates(ν,σ)
    noPositive=true
    v = similar(ν)
    while noPositive
        v = [rand(Normal(d,σ)) for d in ν]
        any(x->x>0,v) ? noPositive=false : nothing
    end
    return v
end

function rand(d::LBA)
    @unpack τ,A,k,ν,σ = d
    b=A+k
    N = length(ν)
    v = sampleDriftRates(ν,σ)
    a = rand(Uniform(0,A),N)
    dt = @. (b-a)/v
    choice,mn = selectWinner(dt)
    rt = τ .+ mn
    return choice,rt
end

function rand(d::LBA,N::Int)
    choice = fill(0,N)
    rt = fill(0.0,N)
    for i in 1:N
        choice[i],rt[i]=rand(d)
    end
    return (choice=choice,rt=rt)
end

logpdf(d::LBA,choice,rt) = log(pdf(d,choice,rt))

function logpdf(d::LBA,data::T) where {T<:NamedTuple}
    return sum(logpdf.(d,data...))
end

function logpdf(dist::LBA,data::Array{<:Tuple,1})
    LL = 0.0
    for d in data
        LL += logpdf(dist,d...)
    end
    return LL
end

function pdf(d::LBA,c,rt)
    @unpack τ,A,k,ν,σ = d
    b=A+k; den = 1.0
    rt < τ ? (return 1e-10) : nothing
    for (i,v) in enumerate(ν)
        if c == i
            den *= dens(d,v,rt)
        else
            den *= (1-cummulative(d,v,rt))
        end
    end
    pneg = pnegative(d)
    den = den/(1-pneg)
    den = max(den,1e-10)
    isnan(den) ? (return 0.0) : (return den)
end

function dens(d::LBA,v,rt)
    @unpack τ,A,k,ν,σ = d
    dt = rt-τ; b=A+k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    dens = (1/A)*(-v*cdf(Normal(0,1),n1) + σ*pdf(Normal(0,1),n1) +
        v*cdf(Normal(0,1),n2) - σ*pdf(Normal(0,1),n2))
    return dens
end

function cummulative(d::LBA,v,rt)
    @unpack τ,A,k,ν,σ = d
    dt = rt-τ; b=A+k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    cm = 1 + ((b-A-dt*v)/A)*cdf(Normal(0,1),n1) -
        ((b-dt*v)/A)*cdf(Normal(0,1),n2) + ((dt*σ)/A)*pdf(Normal(0,1),n1) -
        ((dt*σ)/A)*pdf(Normal(0,1),n2)
    return cm
end

function pnegative(d::LBA)
    @unpack ν,σ=d
    p=1.0
    for v in ν
        p*= cdf(Normal(0,1),-v/σ)
    end
    return p
end

############################################################################################
#                                            DynamicHMC
############################################################################################
struct LBAProb{T}
   data::T
   N::Int
   Nc::Int
end

function (problem::LBAProb)(θ)
    @unpack data=problem
    @unpack v,A,k,tau=θ
    d=LBA(ν=v,A=A,k=k,τ=tau)
    logpdf(d,data)+sum(logpdf.(TruncatedNormal(0,3,0,Inf),v)) +
    logpdf(TruncatedNormal(.8,.4,0,Inf),A)+logpdf(TruncatedNormal(.2,.3,0,Inf),k)+
    logpdf(TruncatedNormal(.4,.1,0,Inf),tau)
end

# Define problem with data and inits.
function sampleDHMC(data,N,Nc,nsamples)
  p = LBAProb(data,N,Nc)
  p((v=fill(.5,Nc),A=.8,k=.2,tau=.4))
  # Write a function to return properly dimensioned transformation.
  problem_transformation(p::LBAProb) =
  as((v=as(Array,asℝ₊,Nc),A=asℝ₊,k=asℝ₊,tau=asℝ₊))
  # Use Flux for the gradient.
  P = TransformedLogDensity(problem_transformation(p), p)
  ∇P = ADgradient(:ForwardDiff, P)
  # FSample from the posterior.
  chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P,nsamples)
  # Undo the transformation to obtain the posterior from the chain.
  posterior = TransformVariables.transform.(Ref(problem_transformation(p)), get_position.(chain));
  return posterior
end
############################################################################################
#                                         Run Code
############################################################################################
Random.seed!(5015)
dist = LBA(ν=[1.0,1.5,2.0],A=.8,k=.2,τ=.4)
N = 10
Nc = 3
data = rand(dist,N)
posterior = sampleDHMC(data,N,Nc,2000)
2

Feel free to ping me (@Tamas_Papp) directly next time. I am very busy this afternoon, but will look into this tonight.

3

Will do. Thank you for your time!

4

Chris I have added a slightly expanded version of your MWE to MCMCBenchmarks/Examples/LBA.

1 Like
5

Thanks Rob. That should be helpful for diagnosing the problem.

6

Can you open an issue at Issues · tpapp/DynamicHMC.jl · GitHub, link this example, and indicate whether you are OK with the package using this code as a test?

7

Thanks for the updated/expanded version. I ran it, and I can’t replicate the error using DynamicHMC v1.0.5; it runs fine. Which version are you using?

8

Thanks for looking into this.

It produces the warning once every 5 or 6 attempts to run it. I’m also using DynamicHMC 1.0.5.

If I enable the Random.seed!(5015) on line 169 it will always fail.

9

I think that adaptation starts in a position that is wacky (discontinuous, or numerically ill-behaved, but I did not debug this). A quick fix is

n = dimension(problem_transformation(p))
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, nsamples; 
                                        q = zeros(n), p = ones(n))

in your sampleDHMC function.

Hope this helps.

1 Like
10

That definitely fixes it! Thanks a lot for your help!

1 Like