Score function Mittag-Leffler pdf

Specific Domains Statistics
1

Good sunday to everyone,

I’ve been trying to compute the score function of the Mittag-Leffler distribution given by the pdf

f(t;\alpha,d)=d t^{\alpha -1}E_{\alpha,\alpha}(-dt^\alpha);\quad d,t>0,\quad \alpha\in(0,1)

where d is known as the rate and \alpha as the fractional parameter. So I’ve been trying to compute

\frac{\partial \log f(t;\alpha,d)}{\partial \alpha}

but I haven’t put a lot of thought onto an efficient way of compute this! I just want to check if something works, so I’ve tried to use ForwardDiff.jl together with MittagLeffler.jl, where they use the numerical methods found in this paper
https://www.researchgate.net/publication/237433354_Computation_of_the_Mittag-Leffler_function_and_its_derivatives
and I’ve been trying to do this

using MittagLeffler, ForwardDiff

d = 127.03238985497151
α = 0.8979932142015109
t = 0.008976795220502681

f(α) = log(d*t^(α-1)*mittleff(α,α,-d*t^α))
ForwardDiff.derivative(f,α) #NaN

I sometimes also get this error with other parameters, I got this error with the parameters above but when I rerun it I just get a NaN don’t know why.

ERROR: MethodError: no method matching kronrod(::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1}}, ::Int64)

Closest candidates are:
  kronrod(::Any, ::Integer, ::Real, ::Real; rtol, quad)
   @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/weightedgauss.jl:90
  kronrod(::Type{T}, ::Integer) where T<:AbstractFloat
   @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:316
  kronrod(::Integer, ::Real, ::Real)
   @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:341
  ...

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:564 [inlined]
  [2] _cachedrule
    @ ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:564 [inlined]
  [3] cachedrule
    @ ~/.julia/packages/QuadGK/OtnWt/src/gausskronrod.jl:569 [inlined]
  [4] do_quadgk(f::MittagLeffler.var"#1#2"{…}, s::Tuple{…}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:7
  [5] (::QuadGK.var"#50#51")(f::Any, s::Any, ::Any)
    @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:253 [inlined]
  [6] handle_infinities(workfunc::QuadGK.var"#50#51"{…}, f::MittagLeffler.var"#1#2"{…}, s::Tuple{…})
    @ QuadGK ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:145
  [7] #quadgk#49
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:252 [inlined]
  [8] quadgk
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:250 [inlined]
  [9] quadgk
    @ ~/.julia/packages/QuadGK/OtnWt/src/adapt.jl:247 [inlined]
 [10] ourquadgk
    @ ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:24 [inlined]
 [11] Kint
    @ ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:27 [inlined]
 [12] Kint
    @ ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:35 [inlined]
 [13] mittleffints(α::ForwardDiff.Dual{…}, β::ForwardDiff.Dual{…}, z::ForwardDiff.Dual{…}, ρ::Float64)
    @ MittagLeffler ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:118
 [14] _mittleff_slow_with_eps(α::ForwardDiff.Dual{…}, β::ForwardDiff.Dual{…}, z::Complex{…}, ρ::Float64)
    @ MittagLeffler ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:197 [inlined]
 [15] _mittleff0(α::ForwardDiff.Dual{…}, β::ForwardDiff.Dual{…}, z::ForwardDiff.Dual{…})
    @ MittagLeffler ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:216
 [16] _mittleff
    @ MittagLeffler ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:200 [inlined]
 [17] mittleff(α::ForwardDiff.Dual{…}, β::ForwardDiff.Dual{…}, z::ForwardDiff.Dual{…})
    @ MittagLeffler ~/.julia/packages/MittagLeffler/8sEkJ/src/mittag_leffler.jl:168
 [18] f(α::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1})
    @ Main ./REPL[5]:1
 [19] derivative(f::typeof(f), x::Float64)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/derivative.jl:14
 [20] top-level scope
    @ REPL[6]:1

I’ve seen this in the repo:
mittleff fails for some arguments. In particular, some of those that evaluate integrals.

So my guess is that it is the package producing the error and not the ForwardDiff, or maybe quadrature rules cannot be forwarddiffed?

Could I have some advice? I’m thinking of this options I have

  1. Try to implement myself the numerical method and hope that automatic differentiation applies.
  2. Try to tackle directly the score function by myself theoretically.

Will automatic differentiation be as fast as me trying to find a quadrature for the score function? Is just using finite differences a bad idea?
There is no literature over score functions of the Mittag-Leffler density!