Incorrect pdf for Dirichlet distribution

[yurivish]

Hello,

I was using Julia to explore statistical distributions and came across what seems to me to be an error in the pdf of the Dirichlet distribution.

julia> using Distributions

julia> beta = Beta(1, 1);

julia> pdf(beta, 0)
1.0

julia> dirichlet = Dirichlet([1, 1]);

julia> pdf(dirichlet, [0, 1])
NaN

My understanding is that the two distributions above are mathematically equivalent, and therefore the results should be the same for both of the pdfs above. The second case uses [0, 1] because the Dirichlet distribution is defined on a simplex.

I believe the behavior of pdf(beta, 0) is correct.

I also noticed a difference in behavior on points outside of the support: the Beta distribution returns zero density, while the Dirichlet function can error or report incredibly large values:

julia> pdf(beta, -0.5)
0.0

julia> pdf(dirichlet, [-0.5, 1.5])
ERROR: DomainError with -0.5:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).

julia> pdf(Dirichlet([1000, 100]), [1, 1])
1.293606521475018e146

[andreasnoack]

You are completely right. The calculation in

Distributions.jl/src/multivariate/dirichlet.jl

Line 177 in 530b9e0

@inbounds s += (a[i] - 1.0) * log(x[i])

should use the xlogy from StatsFuns to get the limit right here and it would also be necessary with a check for being insupport. Furthermore, it looks like the methods should avoid some Float64 hard coding. It would be great if you could prepare a PR.