Rician logpdf inaccurate

[aplavin]

Minimal example: Rician logpdf for moderate parameter values.
Computing Rician(ν=1, σ=0.12) at x=2.

using Distributions: logpdf, Rician
using SpecialFunctions: besselix

ν, σ, x = 1.0, 0.12, 2.0

# 1. Distributions.jl
println("Distributions.jl:  ", logpdf(Rician(ν, σ), x))

# 2. Exact formula with besselix (Wikipedia Rician PDF)
# same (correct) with besselix from Bessels.jl
σ² = σ^2; z = x*ν/σ²
println("Bessel formula:    ", log(x/σ² * besselix(0, z)) - (x-ν)^2/(2σ²))

# 3. BigFloat series (ground truth)
ν_b, σ_b, x_b = big"1.0", big"0.12", big"2.0"
z_b = x_b*ν_b/σ_b^2
I0 = sum(((z_b/2)^k / factorial(big(k)))^2 for k in 0:300)
println("BigFloat series:   ", Float64(log(x_b/σ_b^2) - (x_b^2+ν_b^2)/(2σ_b^2) + log(I0)))

And scipy:

from scipy.stats import rice
print("scipy:", rice.logpdf(2.0, 1.0/0.12, scale=0.12))

Result:

Distributions.jl:  -33.670645484570354   (WRONG, off by 0.497)
Bessel formula:    -33.17342036436976
BigFloat series:   -33.17342036436976
scipy:             -33.17342036436977

[quildtide]

This one is interesting; the implementation here does some transformations on the logpdf of a Noncentral Chi-Square distribution; that function comes from StatsFuns.jl, which gets all of its Noncentral Chi-Square function implementations from R.

I.e. pdf(Rician(ν, σ), x)) is proportional to pdf(NoncentralChisq(2, (ν / σ)²), (x / σ)²)) in theory, with a multiplier of 2y / σ².

The code verbatim right now is:

result = log(2 * abs(x) / σ^2) + logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)

This doesn't look wrong to me in theory, which casts doubt on the accuracy of logpdf(NoncentralChisq) instead. We are looking at λ=69.44 in this case though, which makes me wonder:

using Distributions
using SpecialFunctions: besseli, besselix
ν, σ, x = 1.0, 0.12, 2.0
ncx_pdf(k, λ, x) = 1/2 * exp(-(x+λ)/2) * (x/λ)^(k/4 - 1/2) * besseli(k/2 - 1, sqrt(λ*x))
ncx_logpdf(k, λ, x) = log(1/2) - (x+λ)/2 + (k/4 - 1/2)/log(x/λ) + log(besselix(k/2 - 1, sqrt(λ*x))) + sqrt(λ*x)
println("naive log of pdf: ", log(ncx_pdf(2, (ν / σ)^2, (x / σ)^2)))
println("alt logpdf implementation: ", ncx_logpdf(2, (ν / σ)^2, (x / σ)^2))
println("logpdf(NoncentralChisq) ", logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2))

Yielded:

naive log of pdf: -38.80024179788984
alt logpdf implementation: -38.80024179788984
logpdf(NoncentralChisq) -39.29746691809042

Got the same 0.497 issue here. There's already a motivation to replace the leftover R code; guess we have additional motivation now.

[quildtide]

So in theory, JuliaStats/StatsFuns.jl#202 , an unrelated pull request from last week, fixes this, by coincidence.

But I don't think the Rician distribution needs to rely on Noncentral Chisq either, because there's a lot of extra math since the degrees of freedom is fixed (i.e. besselix(k/2 - 1, x) where k=2 seems silly.

[oscardssmith]

using besseli0 there will be a lot faster if you're willing to pick up a Bessels.jl dependency.

[aplavin]

using besseli0 there will be a lot faster if you're willing to pick up a Bessels.jl dependency.

Yes, I noticed that quite a while ago (#1979), wish Distributions maintainers were more willing to adopt it...

[aplavin]

Wow, perfect timing with that StatsFuns PR :)