add tdist cdfs and quantiles

[oscardssmith]

What should I have for type asserts? They were not added with any specific principles.

[devmotion]

You can have a look at other distributions, e.g., https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/distrs/beta.jl. Here we want to support inputs of type ::Real, e.g., tdistcdf(nu::Real, x::Real) etc. However, similar to the beta.jl file in some cases it might be easiest to first promote both arguments, e.g., tdistcdf(nu::Real, x::Real) = tdistcdf(promote(nu, x)...) and then define tdistcdf(nu::T, x::T) where {T<:Real} since it can make it easier to ensure that constants etc are of the desired type (in this example we would like the result to be of type float(T)).

[devmotion]

This causes undesired promotions. We want e.g. that Float32 inputs yield Float32 results.

[devmotion]

Same here.

[ararslan]

Suggested change

	if x < 0 && ν < (0x1p-26 * x)^2
	return -log(abs(x))*ν + log(ν)*muladd(ν, T(.5), -1.) - logbeta(ν/2, T(.5))
	end
	q = T(0.5)*beta_inc(ν/2, T(.5), ν/muladd(x, x, ν))[1]
	return ifelse(q > 0, 1 - q, q)
	half = one(T) / 2
	if x < 0 && ν < (0x1p-26 * x)^2
	return -log(abs(x)) * ν + log(ν) * muladd(ν, half, -1) - logbeta(ν / 2, half)
	end
	q = half * beta_inc(ν / 2, half, ν / muladd(x, x, ν))[1]
	return ifelse(q > 0, 1 - q, q)

Summary of suggested changes:

• Functional:
  • Change -1. to -1 to avoid promotion to Float64.
  • Use one(T) / 2 in place of T(0.5) since T may not be closed under division, so this ensures that the use of 1/2 matches the type of ν / 2. Alternatively, you could use T(0.5) * ν in place of ν / 2 or 1 // 2 in place of T(0.5). (If the latter, you'd just need to ensure there's a method of beta_inc that can deal with Rationals.)
• Stylistic:
  • The body of the if had an errant 3-space indent.
  • The surrounding code in this file uses spaces around * and / so this should be consistent.

These suggestions apply to the other functions here as well.

[paschermayr]

Thanks a lot for the work! I think in the tdistcdf function, this line

return ifelse(q > 0, 1 - q, q)

should be changed to

return ifelse(x > 0, 1 - q, q)

Else one gets wrong cdfs for negative data, i.e.:

tdistcdf(1., 100.) #0.997
tdistcdf(1., -100.) #0.997

[ararslan]

Suggested change

	return -tdistinvcdf(d, 1-x)
	return -tdistinvcdf(ν, 1 - x)

[paschermayr]

Also, this is only an observation wrt usage of the new cdf, but I think the proposed cdf version does not yet work with Automatic differentiation, as the function beta_inc is not AD compatible. MWE:

using Distributions, DistributionsAD, StatsBase
using ForwardDiff, ReverseDiff
using StatsFuns
function _tdistcdf(ν::S, x::T) where {S<:Real, T<:Real}
    if x < 0 && ν < (0x1p-26 * x)^2
       return -log(abs(x))*ν + log(ν)*StatsFuns.muladd(ν, T(.5), -1.) - StatsFuns.logbeta(ν/2, T(.5))
    end
    q = T(0.5)*StatsFuns.beta_inc(ν/2, T(.5), ν/muladd(x, x, ν))[1]
    return ifelse(x > 0, 1 - q, q)
end
function tdistinvcdf(ν, x::Real)
    if x > .5
        return -tdistinvcdf(d, 1-x)
    end
end

function mytargetfunction(data::AbstractVector)
    function obtaingradient(ν::AbstractVector{R}) where {R<:Real}
        nu = ν[1]
        data_uniform = [_tdistcdf(nu, data[iter]) for iter in eachindex(data)]
        return sum(data_uniform)
    end
end

#working
ν = [3.0]
data = randn(1000)
target = mytargetfunction(data)
target(ν)
#not working
ForwardDiff.gradient(target, ν) #MethodError: no method matching _beta_inc(::ForwardDiff.Dual
ReverseDiff.gradient(target, ν) #MethodError: no method matching _beta_inc(::ReverseDiff.TrackedReal

[ararslan]

Superseded by #149.