# Use FiniteDiff.jl in Turing.jl

**URL:** <https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [May 3, 2024, 9:39pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816 "2024-05-03T21:39:42Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 3, 2024, 9:39pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/1 "2024-05-03T21:39:43Z")

</div>

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:

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

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

```julia
# 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:

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [May 3, 2024, 9:48pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/2 "2024-05-03T21:48:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 3, 2024, 10:07pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/3 "2024-05-03T22:07:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [May 3, 2024, 10:14pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/4 "2024-05-03T22:14:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 4, 2024, 7:37am UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/5 "2024-05-04T07:37:47Z")

</div>

> [@mleprovost](#):
>
> Is there a support for finite difference calculation with [FiniteDiff.jl](https://juliahub.com/ui/Packages/General/FiniteDiff) in [Turing.jl](https://juliahub.com/ui/Packages/General/Turing)? If so, how should I pass it for MAP estimation and HMC/NUTS sampling?

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.

> <https://github.com/tpapp/LogDensityProblemsAD.jl/pull/29>
>
> Fix #26
> 
> This PR removes nearly all of the code in the package extensions to r…ely on \[DifferentiationInterface.jl\](https://github.com/gdalle/DifferentiationInterface.jl) instead.
> 
> \*\*Benefits\*\*
> 
> \- Easier to maintain the AD bindings in one single location
> \- Opens the door to more backends, like Tapir.jl (see https://github.com/withbayes/Tapir.jl/issues/95)
> 
> \*\*Changes (done and planned)\*\*
> 
> \> I did my best to make the changes non-breaking but please double-check
> 
> \- \[x\] Bump Julia compat to 1.10 and package version to 1.10.0
> \- \[x\] Make \`ADgradient\` a concrete struct with 3 fields:
> - \`backend\`: an ADTypes.jl object
> - \`ℓ\`: needs no introduction
> - \`extras\`: stores things like tapes and configs in the DifferentiationInterface.jl paradigm
> \- \[x\] Define translations from \`Symbol\`s to ADTypes.jl in order to preserve every documented constructor (I think, but please double check)
> \- \[x\] Keep special cases in 3 extensions (Enzyme, FiniteDifferences, ForwardDiff) whose constructors need specific types
> \- \[x\] Remove irrelevant tests
> \- \[\] Figure out why Tracker tests are broken
> - "loss is infinite" and failed type inference =\> maybe something to improve in DifferentiationInterfaceTrackerExt?
> \- \[\] Add docstrings

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 4, 2024, 8:38am UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/6 "2024-05-04T08:38:13Z")

</div>

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

> [@mleprovost](#):
>
> ```julia
> ω ~ Normal( 0, 1 ) # standard normal prior 
> x = gammainvccdf(1.0, 1, 0.5*erfc(ω/sqrt(2)) ) 
> 
> ```

is’nt that equivalent to :

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

```julia
x ~ Gamma(1.0, 1) 

```

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

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 6, 2024, 6:39pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/7 "2024-05-06T18:39:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 6, 2024, 6:42pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/8 "2024-05-06T18:42:02Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 6, 2024, 7:03pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/9 "2024-05-06T19:03:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 6, 2024, 7:07pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/10 "2024-05-06T19:07:24Z")

</div>

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](https://turing.ml/v0.22/docs/using-turing/autodiff#switching-ad-modes)

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

> **[GitHub - tpapp/LogDensityProblemsAD.jl: AD backends for LogDensityProblems.jl.](https://github.com/tpapp/LogDensityProblemsAD.jl)**
>
> AD backends for LogDensityProblems.jl. Contribute to tpapp/LogDensityProblemsAD.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 6, 2024, 7:27pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/11 "2024-05-06T19:27:09Z")

</div>

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](https://github.com/tpapp/LogDensityProblems.jl/blob/master/src/LogDensityProblems.jl)

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.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 6, 2024, 7:30pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/12 "2024-05-06T19:30:23Z")

</div>

Absolutely:

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

```

---

<div class="post-metadata">

**Author:** ![devmotion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devmotion/32/9201_2.png) [@devmotion](https://discourse.julialang.org/u/devmotion)\
**Post date:** [May 6, 2024, 8:23pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/13 "2024-05-06T20:23:58Z")

</div>

The page you linked to is outdated, the Turing website is [turinglang.org](http://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](https://turinglang.org/stable/docs/using-turing/autodiff#switching-ad-modes)

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 7, 2024, 5:41am UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/14 "2024-05-07T05:41:42Z")

</div>

My bad, thanks for correcting the link.

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

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 8, 2024, 11:43am UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/15 "2024-05-08T11:43:53Z")

</div>

After correcting StatsFun to StatsFuns, I get

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

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 8, 2024, 1:05pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/16 "2024-05-08T13:05:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 8, 2024, 4:15pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/17 "2024-05-08T16:15:00Z")

</div>

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

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

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

---

<div class="post-metadata">

**Author:** ![janglaubitz](https://avatars.discourse-cdn.com/v4/letter/j/6a8cbe/32.png) [@janglaubitz](https://discourse.julialang.org/u/janglaubitz)\
**Post date:** [May 8, 2024, 9:37pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/18 "2024-05-08T21:37:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [May 8, 2024, 9:55pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/19 "2024-05-08T21:55:54Z")

</div>

> [@janglaubitz](#):
>
> 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.

Glad you like it.

> [@janglaubitz](#):
>
> 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.

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.

> [@janglaubitz](#):
>
> I think resolving AD for gamma\_inv\_cdf in Turing would generally be desired.

Agreed, 100%.

---

<div class="post-metadata">

**Author:** ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)\
**Post date:** [May 8, 2024, 10:07pm UTC](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816/20 "2024-05-08T22:07:07Z")

</div>

> [@lrnv](#):
>
> This is stupid, you should use the gamme prior for faster computation. Note that if seeded correctly, the chain should be _exactly_ the same.

I don’t think this is kind or helpful commentary, nor do I believe it adheres to the [Julia community standards](https://julialang.org/community/standards/#be_respectful_and_inclusive).

[Next page](https://discourse.julialang.org/t/use-finitediff-jl-in-turing-jl/113816.md?page=2)