pdf and rand do not work on Beta with Dual type parameters

[xukai92]

For example the following code fails:

using Distributions
using ForwardDiff: Dual
pdf(Beta(Dual(2.), Dual(2.)), 1)

Error:

ERROR: MethodError: `convert` has no method matching convert(::Type{Float64}, ::ForwardDiff.Dual{0,Float64})
This may have arisen from a call to the constructor Float64(...),
since type constructors fall back to convert methods.
Closest candidates are:
  call{T<:AbstractFloat}(::Type{T<:AbstractFloat}, ::Real, ::RoundingMode{T})
  call{T}(::Type{T}, ::Any)
  convert(::Type{Float64}, ::Int8)
  ...
 in pdf at /Users/kai/.julia/v0.4/Distributions/src/univariates.jl:320

[jmxpearson]

This is a known failure. Beta is defined using _delegate_statsfuns, which defines pdf, cdf, etc. in terms of the functions from StatsFuns.jl. Unfortunately, StatsFuns imports these from Rmath. As a result, they're still only designed for Float64.

It would be great to fix this, but we need StatsFuns to be upgraded. My understanding is that people are working on removing RMath dependency (#294), but that's not finished.

[devmotion]

This is fixed in recent versions of Distributions and StatsFuns:

julia> pdf(Beta(ForwardDiff.Dual(2.), ForwardDiff.Dual(2.)), 1)
Dual{Nothing}(0.0)

julia> pdf(Beta(ForwardDiff.Dual(2.), ForwardDiff.Dual(2.)), 0.5)
Dual{Nothing}(1.5)

[oschulz]

Thanks for this, I ran into that multiple times ... would it be possible to make cdf and quantile of Beta compatible with Dual as well, with what we have? Would be very helpful for variate transformation chains.

[devmotion]

Currently, we use Rmath for cdf and quantile of basically anything that is handled via StatsFuns. There are plans to change this (JuliaStats/StatsFuns.jl#113) and eg for the cdf of the beta distributions we can use SpecialFunctions.beta_inc. However, SpecialFunctions.beta_inc only supports arguments of type T for which float(T) is Float16, Float32 or Float64, so it would still not be possible to pass a ForwardDiff.Dual. I assume one has to add a rule to DiffRules to add support for Dual.

[devmotion]

Related: JuliaMath/SpecialFunctions.jl#337

[oschulz]

Thanks!

[oschulz]

I wonder if we could just allow dual numbers in the beta_inc... functions, at least Dual{Float64} (is there any precision-specific code in there?).

[devmotion]

The implementation is optimized for Float64 (Float32, Float16, and other argument types are only supported by converting them to Float64): https://github.com/JuliaMath/SpecialFunctions.jl/blob/a1df96daaeae3ec663593a7c2d76f4245725ac5e/src/beta_inc.jl#L727-L879

[oschulz]

Ah - but shouldn't Dual{Float64} work too, then?

[oschulz]

Or does the implementation do some low-level bit magic or so?