Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Internals of the incomplete gamma function does not work with ForwardDiff.jl's Dual() type due to enforced floating point precision #465

Open
AkchurinDA opened this issue Feb 19, 2024 · 4 comments

Comments

@AkchurinDA
Copy link

I'm trying to run sensitivity analysis using Fortuna.jl in the example shown here:

using Fortuna

# Define the random vector:
XR      = generaterv("Lognormal", "M", [1.10, 0.10 * 1.10])
XD      = generaterv("Normal",    "M", [1.05, 0.10 * 1.05])
XLApt   = generaterv("Gamma",     "M", [0.22, 0.54 * 0.22])
XLMax   = generaterv("Gumbel",    "M", [1.10, 0.19 * 1.10])
XWMax   = generaterv("Gumbel",    "M", [0.47, 0.35 * 0.47])
X       = [XR, XD, XLApt, XLMax, XWMax]
ρˣ      = Matrix(1.0 * I, 5, 5)

# Define the limit state function:
function g(x, θ)
    R2D = (1 / θ[1]) * max(
        1.4, 
        1.2 + 1.6 * θ[2],
        1.2 + 0.5 * θ[2] + 1.0 * θ[3])

    Q2D = max(
        x[2] + θ[2] * x[3],
        x[2] + θ[2] * x[4],
        x[2] + θ[2] * x[3] + θ[3] * x[5])

    return R2D * x[1] - Q2D
end

# Define the parameters of the limit state function:
θ       = [0.90, 2, 1000] 

# Define and solve the sensitivity problem:
Problem     = SensitivityProblem(X, ρˣ, g, θ)
Solution    = analyze(Problem)

The example fails to run due to the fact that internals of the incomplete gamma function, such as __gamma_inc_inv(), enforce Float64 type and, as the result, errors on ForwardDiff.jl's Dual() type. Here a truncated error message:

ERROR: MethodError: no method matching __gamma_inc_inv(::ForwardDiff.Dual{…}, ::ForwardDiff.Dual{…}, ::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
@AkchurinDA AkchurinDA changed the title Internals of inverse of gamma function does not work with ForwardDiff.jl's Dual() type due to enforced floating point precision Internals of the incomplete gamma function does not work with ForwardDiff.jl's Dual() type due to enforced floating point precision Feb 19, 2024
@stevengj
Copy link
Member

Most special functions are only implemented in a particular precision, so you can't plug in arbitrary numeric types and expect them to work.

In principle, it should be possible to make some such methods work with a type like Dual{Float64} that is based on the same precision, but this would need to be done on a case-by-case basis.

@Deduction42
Copy link

Deduction42 commented Feb 25, 2024

Just to chime in, some of these methods do work with Zygote.jl, I had some provable success using it with __gamma_inc(...) when for the longest time I had to use finite difference. That being said, I can only get gradients with Zygote.jl, because Hessians actually use ForwardDiff.jl (using Zygote.jl on itself has some serious problems); this is a shame because most maximum likelihood problems have very few parameters and are expensive to evaluate, so being able to efficiently and reliably calculate Hessians enables much better optimization algorithms. I can still get by with finite difference methods, but being able to use ForwardDiff.jl on special functions that commonly apply to statistics (particularly anything related to Beta, and Gamma distributions, but I think Normal is covered) opens a lot of doors for very good maximum likelihood estimation solutions.

Maybe Enzyme.jl will solve these issues, but I'm not sure. However, at the time of writing, I can't even get gradients on Gamma cdfs with it.

@AkchurinDA
Copy link
Author

@Deduction42 Thank you for the comment! It is great to know that you can still work it out in Zygote.jl. Would you mind providing a small MWE of how you made it work? I think it will be useful to me and people who stumble upon the issue in the future.

@Deduction42
Copy link

I used this.

using Distributions
using Zygote

f(θ)  = cdf(Gamma(θ[1], θ[2]), 0.5)
∂f(θ) = gradient(f, θ)

display("Calculating gradient...")
∂f([2.0,0.5])

It's actually the first time autodiff has worked for me (last time I tried was three years ago) on this. I tend to run into Bayesian inference and maximum-likelihood estimates on Gamma and Beta distributions from time to time. Zygote has come a long way, and Enzyme is supposed to be even better when it fully arrives.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants