Add Rician distribution

docs/src/univariate.md

@@ -132,6 +132,7 @@ NormalInverseGaussian
 Pareto
 PGeneralizedGaussian
 Rayleigh
+Rician
 Semicircle
 StudentizedRange
 SymTriangularDist

src/Distributions.jl

@@ -142,6 +142,7 @@ export
     PoissonBinomial,
     QQPair,
     Rayleigh,
+    Rician,
     Semicircle,
     Skellam,
     SkewNormal,
@@ -329,7 +330,7 @@ Supported distributions:
     MvNormalKnownCov, MvTDist, NegativeBinomial, NoncentralBeta, NoncentralChisq,
     NoncentralF, NoncentralHypergeometric, NoncentralT, Normal, NormalCanon,
     NormalInverseGaussian, Pareto, PGeneralizedGaussian, Poisson, PoissonBinomial,
-    QQPair, Rayleigh, Skellam, Soliton, StudentizedRange, SymTriangularDist, TDist, TriangularDist,
+    QQPair, Rayleigh, Rician, Skellam, Soliton, StudentizedRange, SymTriangularDist, TDist, TriangularDist,
     Triweight, Truncated, TruncatedNormal, Uniform, UnivariateGMM,
     VonMises, VonMisesFisher, WalleniusNoncentralHypergeometric, Weibull,
     Wishart, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon,

src/univariate/continuous/rician.jl

@@ -0,0 +1,153 @@
+"""
+    Rician(ν, σ)
+
+The *Rician distribution* with parameters `ν` and `σ` has probability density function:
+
+```math
+f(x; \\nu, \\sigma) = \\frac{x}{\\sigma^2} \\exp\\left( \\frac{-(x^2 + \\nu^2)}{2\\sigma^2} \\right) I_0\\left( \\frac{x\\nu}{\\sigma^2} \\right).
+```
+
+If shape and scale parameters `K` and `Ω` are given instead, `ν` and `σ` may be computed from them:
+
+```math
+\\sigma = \\sqrt{\\frac{\\Omega}{2(K + 1)}}, \\quad \\nu = \\sigma\\sqrt{2K}
+```
+
+```julia
+Rician()         # Rician distribution with parameters ν=0 and σ=1
+Rician(ν, σ)     # Rician distribution with parameters ν and σ
+
+params(d)        # Get the parameters, i.e. (ν, σ)
+shape(d)         # Get the shape parameter K = ν²/2σ²
+scale(d)         # Get the scale parameter Ω = ν² + 2σ²
+```
+
+External links:
+
+* [Rician distribution on Wikipedia](https://en.wikipedia.org/wiki/Rice_distribution)
+
+"""
+struct Rician{T<:Real} <: ContinuousUnivariateDistribution
+    ν::T
+    σ::T
+    Rician{T}(ν, σ) where {T} = new{T}(ν, σ)
+end
+
+function Rician(ν::T, σ::T; check_args=true) where {T<:Real}
+    check_args && @check_args(Rician, ν ≥ zero(ν) && σ ≥ zero(σ))
+    return Rician{T}(ν, σ)
+end
+
+Rician() = Rician(0.0, 1.0)
+Rician(ν::Real, σ::Real) = Rician(promote(ν, σ)...)
+Rician(ν::Integer, σ::Integer) = Rician(float(ν), float(σ))
+
+@distr_support Rician 0.0 Inf
+
+#### Conversions
+
+function convert(::Type{Rician{T}}, ν::Real, σ::Real) where T<:Real
+    Rician(T(ν), T(σ))
+end
+
+function convert(::Type{Rician{T}}, d::Rician{S}) where {T <: Real, S <: Real}
+    Rician(T(d.ν), T(d.σ); check_args=false)
+end
+
+#### Parameters
+
+shape(d::Rician) = d.ν^2 / (2 * d.σ^2)
+scale(d::Rician) = d.ν^2 + 2 * d.σ^2
+
+params(d::Rician) = (d.ν, d.σ)
+partype(d::Rician{T}) where {T<:Real} = T
+
+#### Statistics
+
+# helper
+_Lhalf(x) = exp(x/2) * ((1-x) * besseli(zero(x), -x/2) - x * besseli(oneunit(x), -x/2))
+
+mean(d::Rician{T}) where T = d.σ * √T(π/2) * _Lhalf(-d.ν^2/(2 * d.σ^2))
+var(d::Rician) = 2 * d.σ^2 + d.ν^2 - (π * d.σ^2 / 2) * _Lhalf(-d.ν^2/(2 * d.σ^2))^2
+
+quantile(d::Rician, q::Real) = quantile_newton(d, q, mean(d))
+function mode(d::Rician)
+    m = mean(d)
+    _minimize_gss(x -> -pdf(d, x), zero(m), m)
+end
+
+# helper: 1D minimization using Golden-section search
+function _minimize_gss(f, a, b; tol=1e-12)
+    ϕ = (√5 + 1) / 2
+    c = b - (b - a) / ϕ
+    d = a + (b - a) / ϕ
+    while abs(b - a) > tol
+        if f(c) < f(d)
+            b = d
+        else
+            a = c
+        end
+        c = b - (b - a) / ϕ
+        d = a + (b - a) / ϕ
+    end
+    (b + a) / 2
+end
+
+#### PDF/CDF
+
+function logpdf(d::Rician, x::Real)
+    (x ≤ 0 || isinf(x)) && return oftype(float(x), -Inf)
+    σ² = d.σ^2
+    log(x) - log(σ²) - (x^2 + d.ν^2)/(2σ²) + log(besseli(zero(x), x .* (d.ν / σ²)))
+end
+
+function cdf(d::Rician, x::Real)
+    let x = float(x)
+        isnan(x) && return x
+        x ≤ 0.0 && return zero(x)
+        isinf(x) && return one(x)
+        oftype(x, first(quadgk(x -> pdf(d, x), zero(x), x; rtol=1e-8)))
+    end
+end
+
+#### Sampling
+
+function rand(rng::AbstractRNG, d::Rician)
+    x = randn(rng) * d.σ + d.ν
+    y = randn(rng) * d.σ
+    hypot(x, y)
+end
+
+#### Fitting
+
+# implementation based on the Koay inversion technique
+function fit(::Type{<:Rician}, x::AbstractArray{T}; tol=1e-12, maxiters=500) where T
+    μ₁ = mean(x)
+    μ₂ = var(x)
+    r = μ₁ / √μ₂
+    if r < sqrt(π/(4-π))
+        ν = zero(float(T))
+        σ = scale(fit(Rayleigh, float.(x)))
+    else
+        ξ(θ) = 2 + θ^2 - π/8 * exp(-θ^2 / 2) * ((2 + θ^2) * besseli(0, θ^2 / 4) + θ^2 * besseli(1, θ^2 / 4))^2
+        g(θ) = sqrt(ξ(θ) * (1+r^2) - 2)
+        θ = 1.0
+        for j ∈ 1:maxiters
+            θ⁻ = θ
+            θ = g(θ)
+            abs(θ - θ⁻) < tol && break
+        end
+        ξθ = ξ(θ)
+        σ = convert(float(T), sqrt(μ₂ / ξθ))
+        ν = convert(float(T), sqrt(μ₁^2 + (ξθ - 2) * σ^2))
+    end
+    Rician(ν, σ)
+end
+
+# Not implemented:
+#   skewness(d::Rician)
+#   kurtosis(d::Rician)
+#   entropy(d::Rician)
+#   mgf(d::Rician, t::Real)
+#   cf(d::Rician, t::Real)
+

src/univariates.jl

@@ -722,6 +722,7 @@ const continuous_distributions = [
     "logitnormal",    # LogitNormal depends on Normal
     "pareto",
     "rayleigh",
+    "rician",
     "semicircle",
     "skewnormal",
     "studentizedrange",

test/ref/continuous/rician.R

@@ -0,0 +1,23 @@
+
+Rician <- R6Class("Rician",
+    inherit = ContinuousDistribution,
+    public = list(
+        names = c("nu", "sigma"),
+        nu = NA,
+        sigma = NA,
+        initialize = function(n=0, s=1) {
+            self$nu <- n
+            self$sigma <- s
+        },
+        supp = function() { c(0, Inf) },
+        properties = function() {
+            n <- self$nu
+            s <- self$sigma
+            list(scale = n^2 + 2 * s^2,
+                 shape = n^2 / (2 * s^2))
+        },
+        pdf = function(x, log=FALSE) { VGAM::drice(x, self$sigma, self$nu, log=log) },
+        cdf = function(x){ VGAM::price(x, self$sigma, self$nu) },
+        quan = function(v){ VGAM::qrice(v, self$sigma, self$nu) }
+    )
+)

test/ref/continuous_test.lst

@@ -149,6 +149,10 @@ Rayleigh()
 Rayleigh(3.0)
 Rayleigh(8.0)
 
+Rician(1.0, 1.0)
+Rician(5.0, 1.0)
+Rician(10.0, 1.0)
+
 StudentizedRange(2.0, 2.0)
 StudentizedRange(5.0, 10.0)
 StudentizedRange(10.0, 5.0)

test/ref/continuous_test.ref.json

@@ -4050,6 +4050,90 @@
     { "q": 0.90, "x": 17.1677282103148 }
   ]
 },
+{
+  "expr": "Rician(1.0, 1.0)",
+  "dtype": "Rician",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+    "scale": 3,
+    "shape": 0.5
+  },
+  "points": [
+    { "x": 0.58680768479231, "pdf": 0.32597689889069, "logpdf": -1.12092876242478, "cdf": 0.0999999999999981 },
+    { "x": 0.8500705146284, "pdf": 0.42713703077601, "logpdf": -0.850650402069471, "cdf": 0.2 },
+    { "x": 1.06959445774872, "pdf": 0.478590262388661, "logpdf": -0.736910449747733, "cdf": 0.300000000000001 },
+    { "x": 1.27356435184454, "pdf": 0.497262612586662, "logpdf": -0.698636996890676, "cdf": 0.399999999999999 },
+    { "x": 1.47547909178815, "pdf": 0.489049243270422, "logpdf": -0.715292092592867, "cdf": 0.500000000000012 },
+    { "x": 1.6862862997492, "pdf": 0.455970051392373, "logpdf": -0.785328148395677, "cdf": 0.599999999999999 },
+    { "x": 1.91973949265777, "pdf": 0.397730897577084, "logpdf": -0.921979639122226, "cdf": 0.700000000000003 },
+    { "x": 2.2010753309479, "pdf": 0.311617302264047, "logpdf": -1.16597943936394, "cdf": 0.800000000000009 },
+    { "x": 2.6019473978067, "pdf": 0.190248024171792, "logpdf": -1.65942666772506, "cdf": 0.900000000000001 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.58680768479231 },
+    { "q": 0.25, "x": 0.962923340192096 },
+    { "q": 0.50, "x": 1.47547909178815 },
+    { "q": 0.75, "x": 2.05189157517447 },
+    { "q": 0.90, "x": 2.6019473978067 }
+  ]
+},
+{
+  "expr": "Rician(5.0, 1.0)",
+  "dtype": "Rician",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+    "scale": 27,
+    "shape": 12.5
+  },
+  "points": [
+    { "x": 3.83337182711774, "pdf": 0.178067275660391, "logpdf": -1.72559384694809, "cdf": 0.1 },
+    { "x": 4.26739311906375, "pdf": 0.283509558779568, "logpdf": -1.26050943934722, "cdf": 0.199999999999999 },
+    { "x": 4.5808282949603, "pdf": 0.351696201209667, "logpdf": -1.04498754078411, "cdf": 0.299999999999981 },
+    { "x": 4.84891183432576, "pdf": 0.390461020883228, "logpdf": -0.940427133165443, "cdf": 0.400000000000014 },
+    { "x": 5.0996760375676, "pdf": 0.402913231155499, "logpdf": -0.909034047523853, "cdf": 0.499999999999981 },
+    { "x": 5.35060611077152, "pdf": 0.389944217738904, "logpdf": -0.941751581527125, "cdf": 0.599999999999959 },
+    { "x": 5.61923787901142, "pdf": 0.350723907317018, "logpdf": -1.04775595388, "cdf": 0.699999999999992 },
+    { "x": 5.93381683616078, "pdf": 0.282226926437332, "logpdf": -1.26504382796596, "cdf": 0.799999999999998 },
+    { "x": 6.37038420267928, "pdf": 0.17678589781236, "logpdf": -1.73281589546462, "cdf": 0.899999999999982 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 3.83337182711774 },
+    { "q": 0.25, "x": 4.43248557648109 },
+    { "q": 0.50, "x": 5.0996760375676 },
+    { "q": 0.75, "x": 5.76805276699816 },
+    { "q": 0.90, "x": 6.37038420267928 }
+  ]
+},
+{
+  "expr": "Rician(10.0, 1.0)",
+  "dtype": "Rician",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+    "scale": 102,
+    "shape": 50
+  },
+  "points": [
+    { "x": 8.77189934889994, "pdf": 0.17602367421233, "logpdf": -1.73713678041993, "cdf": 0.0999999999999889 },
+    { "x": 9.21055856310058, "pdf": 0.280747347661651, "logpdf": -1.27030013273982, "cdf": 0.200000000000031 },
+    { "x": 9.52691154521682, "pdf": 0.348625118414454, "logpdf": -1.05375809337315, "cdf": 0.299999999999978 },
+    { "x": 9.79725334900319, "pdf": 0.387340667844654, "logpdf": -0.948450694502043, "cdf": 0.40000000000009 },
+    { "x": 10.0499586252229, "pdf": 0.399938412039171, "logpdf": -0.91644471363081, "cdf": 0.500000000000001 },
+    { "x": 10.3026851248306, "pdf": 0.387275589795634, "logpdf": -0.948618721053336, "cdf": 0.599999999999961 },
+    { "x": 10.5730967493158, "pdf": 0.348503584620269, "logpdf": -1.05410676298003, "cdf": 0.700000000000032 },
+    { "x": 10.8895938007392, "pdf": 0.280589475189091, "logpdf": -1.27086262025283, "cdf": 0.799999999999942 },
+    { "x": 11.3285661035151, "pdf": 0.175871273888067, "logpdf": -1.73800294990951, "cdf": 0.900000000000024 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 8.77189934889994 },
+    { "q": 0.25, "x": 9.37722801591067 },
+    { "q": 0.50, "x": 10.0499586252229 },
+    { "q": 0.75, "x": 10.7228400133078 },
+    { "q": 0.90, "x": 11.3285661035151 }
+  ]
+},
 {
   "expr": "StudentizedRange(2.0, 2.0)",
   "dtype": "StudentizedRange",

test/ref/rdistributions.R

@@ -68,6 +68,7 @@ source("continuous/normal.R")
 source("continuous/normalinversegaussian.R")
 source("continuous/pareto.R")
 source("continuous/rayleigh.R")
+source("continuous/rician.R")
 source("continuous/studentizedrange.R")
 source("continuous/tdist.R")
 source("continuous/triangulardist.R")

test/rician.jl

@@ -0,0 +1,40 @@
+@testset "Rician" begin
+
+    d1 = Rician(0.0, 10.0)
+    @test d1 isa Rician{Float64}
+    @test params(d1) == (0.0, 10.0)
+    @test shape(d1) == 0.0
+    @test scale(d1) == 200.0
+    @test partype(d1) === Float64
+    @test eltype(d1) === Float64
+    @test rand(d1) isa Float64
+
+    d2 = Rayleigh(10.0)
+    @test mean(d1) ≈ mean(d2)
+    @test var(d1) ≈ var(d2)
+    @test mode(d1) ≈ mode(d2)
+    @test median(d1) ≈ median(d2)
+    @test quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) ≈ quantile.(d2, [0.25, 0.45, 0.60, 0.80, 0.90])
+    @test pdf.(d1, 0.0:0.1:1.0) ≈ pdf.(d2, 0.0:0.1:1.0)
+    @test cdf.(d1, 0.0:0.1:1.0) ≈ cdf.(d2, 0.0:0.1:1.0)
+
+    x = rand(Rician(5.0, 5.0), 100000)
+    d1 = fit(Rician, x)
+    @test d1 isa Rician{Float64}
+    @test params(d1)[1] ≈ 5.0 atol=1e-1
+    @test params(d1)[2] ≈ 5.0 atol=1e-1
+
+    d1 = Rician(10.0f0, 10.0f0)
+    @test d1 isa Rician{Float32}
+    @test params(d1) == (10.0f0, 10.0f0)
+    @test shape(d1) == 0.5f0
+    @test scale(d1) == 300.0f0
+    @test partype(d1) === Float32
+    @test eltype(d1) === Float64
+    @test rand(d1) isa Float64
+
+    d1 = Rician()
+    @test d1 isa Rician{Float64}
+    @test params(d1) == (0.0, 1.0)
+
+end

test/runtests.jl

@@ -62,6 +62,7 @@ const tests = [
     "chi",
     "gumbel",
     "pdfnorm",
+    "rician",
 ]
 
 printstyled("Running tests:\n", color=:blue)