[Largely edited version]
I was looking at this issue about the derivative of the incomplete gamma function w.r.t the parameter a.
Looks like from this page of wolfram alpha docs that the right formula is :
![]()
which depends on the regularized hypergeometric function. After a certain time of translation, namely using this page to match the wolfram alpha defintion to the one from the HypegeometricFunctions.jl package (also, the formula is about non-normalized upper incomplete gamma function, while gamma_inc gives the normalized version), I have found a working implementation :
using SpecialFunctions, HypergeometricFunctions
function Γ(a,x)
return gamma_inc(a,x)[2]
end
function ∂Γ_∂a(a,x)
g = gamma(a)
dg = digamma(a)
G = Γ(a,x)
lx = log(x)
r = pFq([a,a],[a+1,a+1],-x)
return a^(-2) * x^a * r/g + (1 - G)*(dg - lx)
end
function emp(a,x)
ϵ = 1e-5
return (Γ(a+ϵ,x) - Γ(a-ϵ,x))/2ϵ
end
a = 1.9
x = 1.8
Γ(a,x)
(∂Γ_∂a(a,x), emp(a,x))
#a simple test:
a = 0.1:0.1:5
x = 0.1:0.1:5
diff(a,x) = abs(∂Γ_∂a(a,x) - emp(a,x))
maximum(diff.(a,x')) # 2e-10 on my machine.
How could I:
- Check in a better way that this is yielding the right derivative ?
- Tell ForwardDiff to use this instead of erroring out on
ForwardDiff.gradient(a -> gamma_inc(a,1)[2],1.2).