Use FiniteDiff.jl in Turing.jl

Specific Domains Probabilistic programming
1

Hello,

I have a user-defined density which does not support automatic differentiation (the function gamma_inv_cdf in SpecialFunctions.jl does not support AD).

Is there a support for finite difference calculation with FiniteDiff.jl in Turing.jl? If so, how should I pass it for MAP estimation and HMC/NUTS sampling?

As an alternative, can I provide analytical calcuation of the gradient of the log-density?

MWI:

using Turing
using SpecialFunctions
using StatsFun

# Define the model for the original toy model.  
@model function model_toyProblem_tranformed(y)
    ω ~ Normal( 0, 1 ) # standard normal prior 
    x = gammainvccdf(1.0, 1, 0.5*erfc(ω/sqrt(2)) ) 

    y ~ Normal( x, 1.0) # pushforward likelihood 
end

model_func = model_toyProblem_tranformed(0.0)

map_estimate = optimize( model_func, MAP()); # MAP estimate

Output:

MethodError: no method matching __gamma_inc_inv(::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}, ::Bool)

Closest candidates are:
  __gamma_inc_inv(::Float64, ::Float64, ::Bool)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:1012
  __gamma_inc_inv(::T, ::T, ::Bool) where T<:Union{Float16, Float32}
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:1089


Stacktrace:
  [1] _gamma_inc_inv
    @ ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:1009 [inlined]
  [2] gamma_inc_inv
    @ ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:991 [inlined]
  [3] gammainvccdf
    @ ~/.julia/packages/StatsFuns/mrf0e/src/distrs/gamma.jl:108 [inlined]
  [4] model_toyProblem_tranformed(__model__::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}, Vector{Set{DynamicPPL.Selector}}}}, ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}, __context__::Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}, y::Float64)
    @ Main ./In[15]:4
  [5] _evaluate!!
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/model.jl:963 [inlined]
  [6] evaluate_threadunsafe!!
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/model.jl:936 [inlined]
  [7] evaluate!!
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/model.jl:889 [inlined]
  [8] (::LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}})(z::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}})
    @ Turing.Optimisation ~/.julia/packages/Turing/lkUBK/src/optimisation/Optimisation.jl:111
  [9] logdensity
    @ ~/.julia/packages/Turing/lkUBK/src/optimisation/Optimisation.jl:115 [inlined]
 [10] Fix1
    @ ./operators.jl:1118 [inlined]
 [11] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24 [inlined]
 [12] vector_mode_gradient!
    @ ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:96 [inlined]
 [13] gradient!(result::DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, f::Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}}}, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:37
 [14] gradient!
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:35 [inlined]
 [15] logdensity_and_gradient(fℓ::LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}}, ForwardDiff.Chunk{1}, ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}}, x::Vector{Float64})
    @ LogDensityProblemsADForwardDiffExt ~/.julia/packages/LogDensityProblemsAD/rBlLq/ext/LogDensityProblemsADForwardDiffExt.jl:118
 [16] (::LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}})(F::Float64, G::Vector{Float64}, z::Vector{Float64})
    @ Turing.Optimisation ~/.julia/packages/Turing/lkUBK/src/optimisation/Optimisation.jl:122
 [17] (::NLSolversBase.var"#69#70"{NLSolversBase.InplaceObjective{Nothing, LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}}, Nothing, Nothing, Nothing}, Float64})(G::Vector{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/incomplete.jl:54
 [18] value_gradient!!(obj::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:82
 [19] initial_state(method::LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}, options::Optim.Options{Float64, Nothing}, d::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64})
    @ Optim ~/.julia/packages/Optim/DV9Sq/src/multivariate/solvers/first_order/l_bfgs.jl:164
 [20] optimize(d::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}, options::Optim.Options{Float64, Nothing}, state::Optim.LBFGSState{Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Float64, Vector{Float64}})
    @ Optim ~/.julia/packages/Optim/DV9Sq/src/multivariate/optimize/optimize.jl:36 [inlined]
 [21] optimize(f::NLSolversBase.InplaceObjective{Nothing, LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}}, Nothing, Nothing, Nothing}, initial_x::Vector{Float64}, method::LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}, options::Optim.Options{Float64, Nothing}; inplace::Bool, autodiff::Symbol)
    @ Optim ~/.julia/packages/Optim/DV9Sq/src/multivariate/optimize/interface.jl:143
 [22] optimize
    @ ~/.julia/packages/Optim/DV9Sq/src/multivariate/optimize/interface.jl:139 [inlined]
 [23] _optimize(::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, ::LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}}, ::Vector{Float64}, ::LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}, ::Optim.Options{Float64, Nothing}; kwargs::@Kwargs{})
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:235
 [24] _optimize(::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, ::LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, Turing.Optimisation.OptimizationContext{DynamicPPL.DefaultContext}}, ::Vector{Float64}, ::LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Optim.var"#19#21"}, ::Optim.Options{Float64, Nothing})
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:219
 [25] _map_optimize(::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, ::Vector{Float64}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:211
 [26] _map_optimize(::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, ::Vector{Float64}, ::Vararg{Any})
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:209
 [27] optimize(model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, ::MAP, options::Optim.Options{Float64, Nothing}; kwargs::@Kwargs{})
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:186
 [28] optimize
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:181 [inlined]
 [29] optimize(model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, ::MAP)
    @ TuringOptimExt ~/.julia/packages/Turing/lkUBK/ext/TuringOptimExt.jl:181
 [30] top-level scope
    @ In[18]:1

For NUTS sampling (note the double S)

# Prepare arguments. 
nr_samples = 10^4 # number of samples

# Generate a chain containing the samples of x and θ
chn_NUTS = sample(model_func, 
    NUTS(), 
    nr_samples)

Output:

MethodError: no method matching __gamma_inc_inv(::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}, ::Bool)

Closest candidates are:
  __gamma_inc_inv(::Float64, ::Float64, ::Bool)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:1012
  __gamma_inc_inv(::T, ::T, ::Bool) where T<:Union{Float16, Float32}
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:1089


Stacktrace:
  [1] _gamma_inc_inv
    @ ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:1009 [inlined]
  [2] gamma_inc_inv
    @ ~/.julia/packages/SpecialFunctions/QH8rV/src/gamma_inc.jl:991 [inlined]
  [3] gammainvccdf
    @ ~/.julia/packages/StatsFuns/mrf0e/src/distrs/gamma.jl:108 [inlined]
  [4] model_toyProblem_tranformed(__model__::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}, Vector{Set{DynamicPPL.Selector}}}}, ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}, y::Float64)
    @ Main ./In[15]:4
  [5] _evaluate!!
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/model.jl:963 [inlined]
  [6] evaluate_threadunsafe!!
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/model.jl:936 [inlined]
  [7] evaluate!!
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/model.jl:889 [inlined]
  [8] logdensity(f::LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}, θ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/rXg4T/src/logdensityfunction.jl:94
  [9] Fix1
    @ ./operators.jl:1118 [inlined]
 [10] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24 [inlined]
 [11] vector_mode_gradient!
    @ ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:96 [inlined]
 [12] gradient!(result::DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, f::Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}}, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:37
 [13] gradient!(result::DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, f::Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}}, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:35 [inlined]
 [14] logdensity_and_gradient
    @ ~/.julia/packages/LogDensityProblemsAD/rBlLq/ext/LogDensityProblemsADForwardDiffExt.jl:118 [inlined]
 [15] ∂logπ∂θ
    @ ~/.julia/packages/Turing/lkUBK/src/mcmc/hmc.jl:159 [inlined]
 [16] ∂H∂θ
    @ ~/.julia/packages/AdvancedHMC/AlvV4/src/hamiltonian.jl:38 [inlined]
 [17] phasepoint(h::AdvancedHMC.Hamiltonian{AdvancedHMC.DiagEuclideanMetric{Float64, Vector{Float64}}, AdvancedHMC.GaussianKinetic, Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}, ForwardDiff.Chunk{1}, ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}}}, Turing.Inference.var"#∂logπ∂θ#31"{LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}, ForwardDiff.Chunk{1}, ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}}}}, θ::Vector{Float64}, r::Vector{Float64})
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/AlvV4/src/hamiltonian.jl:74
 [18] phasepoint(rng::Random.TaskLocalRNG, θ::Vector{Float64}, h::AdvancedHMC.Hamiltonian{AdvancedHMC.DiagEuclideanMetric{Float64, Vector{Float64}}, AdvancedHMC.GaussianKinetic, Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}, ForwardDiff.Chunk{1}, ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}}}, Turing.Inference.var"#∂logπ∂θ#31"{LogDensityProblemsADForwardDiffExt.ForwardDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random.TaskLocalRNG}}, ForwardDiff.Chunk{1}, ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 1}}}}}})
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/AlvV4/src/hamiltonian.jl:155
 [19] initialstep(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, vi_original::DynamicPPL.TypedVarInfo{@NamedTuple{ω::DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ω, Setfield.IdentityLens}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:ω, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}, Float64}; initial_params::Nothing, nadapts::Int64, kwargs::@Kwargs{})
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/hmc.jl:163
 [20] step(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}; initial_params::Nothing, kwargs::@Kwargs{nadapts::Int64})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/rXg4T/src/sampler.jl:116
 [21] step
    @ ~/.julia/packages/DynamicPPL/rXg4T/src/sampler.jl:99 [inlined]
 [22] macro expansion
    @ ~/.julia/packages/AbstractMCMC/YrmkI/src/sample.jl:130 [inlined]
 [23] macro expansion
    @ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
 [24] (::AbstractMCMC.var"#22#23"{Bool, String, Nothing, Int64, Int64, Nothing, @Kwargs{nadapts::Int64}, Random.TaskLocalRNG, DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64})()
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/YrmkI/src/logging.jl:12
 [25] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging ./logging.jl:515
 [26] with_logger
    @ Base.CoreLogging ./logging.jl:627 [inlined]
 [27] with_progresslogger(f::Function, _module::Module, logger::Logging.ConsoleLogger)
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/YrmkI/src/logging.jl:36
 [28] macro expansion
    @ ~/.julia/packages/AbstractMCMC/YrmkI/src/logging.jl:11 [inlined]
 [29] mcmcsample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, initial_state::Nothing, kwargs::@Kwargs{nadapts::Int64})
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/YrmkI/src/sample.jl:120
 [30] sample(rng::Random.TaskLocalRNG, model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; chain_type::Type, resume_from::Nothing, initial_state::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::@Kwargs{})
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/hmc.jl:117
 [31] sample
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/hmc.jl:86 [inlined]
 [32] #sample#3
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/Inference.jl:219 [inlined]
 [33] sample
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/Inference.jl:212 [inlined]
 [34] #sample#2
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/Inference.jl:209 [inlined]
 [35] sample(model::DynamicPPL.Model{typeof(model_toyProblem_tranformed), (:y,), (), (), Tuple{Float64}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{AutoForwardDiff{0, Nothing}, (), AdvancedHMC.DiagEuclideanMetric}, N::Int64)
    @ Turing.Inference ~/.julia/packages/Turing/lkUBK/src/mcmc/Inference.jl:203
 [36] top-level scope
    @ In[19]:7

Thank you for your help,

2

Try to use a nativr implementation of the Gamma functions such as the one in Bessels.jl i think it is not exported.

3

I don’t see support for AD in Bessels.jl. The routines are only defined for Float64 or Float 32 numbers. The code does not contain routine for the inverse of the incomplete Gamma function.

4

The Master branch is under heavy work for AD. I don’t know about the inverse gamma though.

5

I think all gradient calculations for Turing.jl are based on LogDensityProblemsAD.jl, which supports FiniteDifferences.jl but not FiniteDiff.jl at the moment.

I have submitted a PR enabling support for every AD backend out there, but it’s not obvious due to Julia compatibility requirements, and it might have to wait until the next LTS.

1 Like
6

Maybe it is htat your MWE is too simple, but:

is’nt that equivalent to :

u ~ Uniform(0,1)
x = gammainvccdf(1.0, 1, 0.5*u)

where u are the quantiles of your Gaussian (you can then move from u to omega after the MCMC finished via ω = quantile(Normal(0,1), u) ? Maybe there is a scaling factor that is not exactly right, but erfc(x/sqrt(2)) looks a lot like the normal cdf.

Then, isnt that equivalent again to

x ~ Gamma(1.0, 1)

and then simply compute u after as twice the cdf of the realizations ?

7

Good point, but the real setting is actually more complicated, as we are using the generalized Gamma distribution. Thus, we cannot apply your suggestion.

8

We are fine working with FiniteDifferences.jl. Do you know how to say to Turing to use FiniteDifferences.jl for the differentiation?

9

You are sure ? Would you mind sharing the complete problem ?

1 Like
10

I’m not sure, cause the Turing docs says it supports only 4 backends

https://turing.ml/v0.22/docs/using-turing/autodiff#switching-ad-modes

but the underlying package supports 5, the last one being FiniteDifferences

1 Like
11

Is there a way to just pass a user-defined object of type LogDensityOrder{1} to Turing.jl in LogDensityProblems.jl/src/LogDensityProblems.jl at master · tpapp/LogDensityProblems.jl · GitHub

That will bypass the issue of dealing with the choice of differentiation backend in Turing.jl, and even allow providing an analytical version of the log density and score.

12

Absolutely:

using Turing
using StatsFuns
using Optim
using SpecialFunctions

# Fix model parameters
y = 0.5 # observational data 
σ² = 10^(-2) # noise variance 

## Selecion of hyper-prior parameters
# power parameter
r = -1.0
# shape parameter 
β = 1.0017
# rate parameters 
ϑ = 1.2308e-4

# Define the first component of the inverse transport map, T1
function priorNormalizing_transportMap_inv_T1( ω; r, β, ϑ ) 
    R = 35 # value above which we extend T1 linearly 

    if ω <= R         
        T1 = gammainvccdf( β, 1, 0.5*erfc(ω/sqrt(2)) ) 
        T1 = ϑ * T1.^( 1/abs(r) ) # compute the inverse CDF
    else 
        # Value of T1 at ω = R          
        T1_R = gammainvccdf( β, 1, 0.5*erfc(R/sqrt(2)) )
        T1_R = ϑ * T1_R.^( 1/abs(r) ) # compute the inverse CDF 
        
        # Compute the value of the derivative of T1 at ω = R using a FD approximation 
        δ = .1 
        T1_aux = gammainvccdf( β, 1, 0.5*erfc((R-δ)/sqrt(2)) ) 
        T1_aux = ϑ * T1_aux.^( 1/abs(r) ) # compute the inverse CDF
        T1_deriv = ( T1_R - T1_aux ) / δ

        # compute the value of the linear extension 
        T1 = T1_R + T1_deriv.*( ω - R )
    end

    return T1
end


# Define the second component of the inverse transport map, T2
function priorNormalizing_transportMap_inv_T2( v, ω; r, β, ϑ ) 
    θ = priorNormalizing_transportMap_inv_T1( ω; r, β, ϑ ) # compute hyper-parameter value θ 
    T2 = v .* sqrt.(θ)
    return T2
end

# Define the model for the original toy model.  
@model function model_toyProblem_tranformed( y, σ², r, β, ϑ )
    v ~ Normal( 0, 1 ) # standard normal prior 
    ω ~ Normal( 0, 1 ) # standard normal prior 
    x = priorNormalizing_transportMap_inv_T2( v, ω; r, β, ϑ )
    y ~ Normal( x, sqrt(σ²) ) # pushforward likelihood 
end

# Condition the model on the given data and parameters 
model_func = model_toyProblem_tranformed( y, σ², r, β, ϑ )

ml_estimate = optimize( model_func, MLE() ); # ML estimate
13

The page you linked to is outdated, the Turing website is turinglang.org (long story short, nobody has access to turing.ml anymore and someone else aquired it; requests for it to be taken down have been unsuccessful AFAIK).

The correct docs are Automatic Differentiation

2 Likes
14

My bad, thanks for correcting the link.

As for the backends, any reason why LogDensityProblemsAD supports FiniteDifferences but Turing doesnt?

15

After correcting StatsFun to StatsFuns, I get

ml_estimate = optimize( model_func, MLE() ); # ML estimate 
ERROR: UndefVarError: `optimize` not defined
Stacktrace:
 [1] top-level scope
   @ Untitled-1:1

Are other packages needed?

Also, note that

julia> ccdf(Normal(big(0),big(1)),big(35))
1.12491070647240624397924276753507491800121497812382401152901023921682165333983e-268

julia>

So your else statement will rarelly apply… Maybe we can simply avoid implementing it, but since I cannot run your code yet i do not undertand why you needed that.

1 Like
16

Hello,

I forgot to add using Optim
I have corrected the dependencies in the code above.
We can ignore the if-else statement for now, the issue comes from differentiating gammainvccdf not the if-else statement.

17

As i thought, using directly the Gamma distribution instead of its quantile, or even worse a gaussian r.v. works perfectly fine:

using Turing
using SpecialFunctions
using StatsFuns
using Optim

y = 0.5 # observational data 
σ² = 10^(-2) # noise variance 
r = -1.0 # power parameter
β = 1.0017 # shape parameter 
ϑ = 1.2308e-4 # rate parameters 

# Define the model for the original toy model.  
@model function model_toyProblem_tranformed( y, σ², r, β, ϑ )
    v ~ Normal( 0, 1 ) # standard normal prior 
    # ω ~ Normal( 0, 1 ) # standard normal prior 
    # q ~ Uniform(0,1) # standard uniform prior? q = 0.5 * erfc(ω/sqrt(2))
    Γ ~ Gamma(β,1) # direct gamma prior, Γ = quantile(Gamma(β,1), 1 - q) = gamma_inv_inc(...) you had before. 
    x  = v * sqrt(ϑ * Γ^(1/abs(r)))
    y ~ Normal( x, sqrt(σ²) ) # pushforward likelihood 
end
model_func = model_toyProblem_tranformed( y, σ², r, β, ϑ )
ml_estimate = optimize( model_func, MLE() ); # ML estimate

simply outputs :

ModeResult with maximized lp of 1.38
[5.984977813674527, 56.70579342723415]

which I do not understand / tried to interpret, but at least no more error of differentiability of gamma_inc_inv. Few comments:

  • Another solution would have been to define the generalized gamma distribution you want and provide its density directly (which is on wikipedia, but not already inside Distributions.jl weirdly).

  • Note that the map between \omega, q, and \Gamma is bijective so if you want to interpret your results on the gaussian scale you still can, just run the optimization on the natural scale

  • Another good reason to use this version is that it is probably waaaay faster than AD of gamma_inc_inv :wink:

1 Like
18

Hi Irnv,

Thanks for the suggestion. We agree that reformulating the model in this way allows us to avoid the AD problem with gamma_inc_inv for finding ML and MAP estimates.

However, we also want to sample from the model (using NUTS), and there, we strongly prefer staying with v, \omega \sim Normal(0,1) as this corresponds to a standard normal prior.

I think resolving AD for gamma_inv_cdf in Turing would generally be desired.

Thanks,
Jan

1 Like
19

Glad you like it.

This is stupid, you should use the gamme prior for faster computation. Note that if seeded correctly, the chain should be exactly the same. As I said, the clear bijection between the two gives you a normal prior for later interpretation and analysis if you really need it.

Agreed, 100%.

20

I don’t think this is kind or helpful commentary, nor do I believe it adheres to the Julia community standards.

1 Like