From bf5293491257e7a38b6a44af28b1d140b74d9994 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 01:22:16 +0800 Subject: [PATCH 01/24] Add Rician distribution --- docs/src/univariate.md | 1 + src/Distributions.jl | 3 +- src/univariate/continuous/rician.jl | 152 ++++++++++++++++++++++++++++ src/univariates.jl | 1 + test/ref/continuous/rician.R | 23 +++++ test/ref/continuous_test.lst | 4 + test/ref/continuous_test.ref.json | 84 +++++++++++++++ test/ref/rdistributions.R | 1 + test/rician.jl | 47 +++++++++ test/runtests.jl | 1 + 10 files changed, 316 insertions(+), 1 deletion(-) create mode 100644 src/univariate/continuous/rician.jl create mode 100644 test/ref/continuous/rician.R create mode 100644 test/rician.jl diff --git a/docs/src/univariate.md b/docs/src/univariate.md index 0418a38fa..b2bdf6290 100644 --- a/docs/src/univariate.md +++ b/docs/src/univariate.md @@ -132,6 +132,7 @@ NormalInverseGaussian Pareto PGeneralizedGaussian Rayleigh +Rician Semicircle StudentizedRange SymTriangularDist diff --git a/src/Distributions.jl b/src/Distributions.jl index 72005e7f8..52b1b84a3 100644 --- a/src/Distributions.jl +++ b/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, diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl new file mode 100644 index 000000000..aa1a80458 --- /dev/null +++ b/src/univariate/continuous/rician.jl @@ -0,0 +1,152 @@ +""" + 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). +``` + +```julia +Rician() # Rician distribution with parameters ν=0 and σ=1 +Rician(ν, σ) # Rician distribution with parameters ν and σ +Rician(; K, Ω) # Rician distribution with shape K and scale Ω + +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(ν::Real, σ::Real) = Rician(promote(ν, σ)...) +Rician(ν::Integer, σ::Integer) = Rician(float(ν), float(σ)) + +function Rician(; K=0.0, Ω=2.0) + σ = sqrt(Ω / (K + 1) / 2) + ν = √2 * σ * sqrt(K) + Rician(ν, σ) +end + +@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.σ) +@inline partype(d::Rician{T}) where {T<:Real} = T + +#### Statistics + +# helper +L½(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) * L½(-d.ν^2/(2 * d.σ^2)) +var(d::Rician) = 2 * d.σ^2 + d.ν^2 - (π * d.σ^2 / 2) * L½(-d.ν^2/(2 * d.σ^2))^2 + +quantile(d::Rician, q::Real) = quantile_newton(d, q, mean(d)) +mode(d::Rician{T}) where T = _minimize_gss(x -> -pdf(d, x), zero(T), mean(d)) + +# 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) + isnan(x) && return x + x ≤ 0.0 && return -convert(eltype(x), Inf) + isinf(x) && return -convert(eltype(x), Inf) + σ² = d.σ^2 + log.(x) - log(σ²) - (x.^2 + d.ν^2)/(2σ²) + log.(besseli.(zero(eltype(x)), x .* (d.ν / σ²))) +end + +pdf(d::Rician, x::Real) = exp.(logpdf(d, x)) + +function cdf(d::Rician, x::Real) + isnan(x) && return x + x ≤ 0.0 && return zero(eltype(x)) + isinf(x) && return one(eltype(x)) + first(quadgk(x -> pdf(d, x), zero(eltype(x)), x; rtol=1e-8)) +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(T) + σ = scale(fit(Rayleigh, x)) + else + ξ(θ) = 2 + θ^2 - T(π)/8 * exp(-θ^2 / 2) * ((2 + θ^2) * besseli(0, θ^2 / 4) + θ^2 * besseli(1, θ^2 / 4))^2 + g(θ) = sqrt(ξ(θ) * (1+r^2) - 2) + θ = one(T) + for j ∈ 1:maxiters + θ⁻ = θ + θ = g(θ) + abs(θ - θ⁻) < tol && break + end + ξθ = ξ(θ) + σ = sqrt(μ₂ / ξθ) + ν = 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) + diff --git a/src/univariates.jl b/src/univariates.jl index f0b992087..0b232b6dd 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -722,6 +722,7 @@ const continuous_distributions = [ "logitnormal", # LogitNormal depends on Normal "pareto", "rayleigh", + "rician", "semicircle", "skewnormal", "studentizedrange", diff --git a/test/ref/continuous/rician.R b/test/ref/continuous/rician.R new file mode 100644 index 000000000..3214f360b --- /dev/null +++ b/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) } + ) +) diff --git a/test/ref/continuous_test.lst b/test/ref/continuous_test.lst index eaef0e4e0..381a1c466 100644 --- a/test/ref/continuous_test.lst +++ b/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) diff --git a/test/ref/continuous_test.ref.json b/test/ref/continuous_test.ref.json index dd8807a09..bb7ce9a89 100644 --- a/test/ref/continuous_test.ref.json +++ b/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", diff --git a/test/ref/rdistributions.R b/test/ref/rdistributions.R index 699f1f718..c8b3b1031 100644 --- a/test/ref/rdistributions.R +++ b/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") diff --git a/test/rician.jl b/test/rician.jl new file mode 100644 index 000000000..c75318432 --- /dev/null +++ b/test/rician.jl @@ -0,0 +1,47 @@ +@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) + + d1 = Rician(K=0.5, Ω=300.0) + @test d1 isa Rician{Float64} + @test params(d1)[1] ≈ 10.0 + @test params(d1)[2] ≈ 10.0 + @test shape(d1) ≈ 0.5 + @test scale(d1) ≈ 300.0 + +end diff --git a/test/runtests.jl b/test/runtests.jl index e2a83166a..a51b8c04c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,6 +62,7 @@ const tests = [ "chi", "gumbel", "pdfnorm", + "rician", ] printstyled("Running tests:\n", color=:blue) From 9ddff52f1c7c2567d5c23e5aeae0380733103574 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 19:47:28 +0800 Subject: [PATCH 02/24] Removed Rician Integer constructor based on review Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index aa1a80458..15b5affe9 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -34,7 +34,6 @@ function Rician(ν::T, σ::T; check_args=true) where {T<:Real} end Rician(ν::Real, σ::Real) = Rician(promote(ν, σ)...) -Rician(ν::Integer, σ::Integer) = Rician(float(ν), float(σ)) function Rician(; K=0.0, Ω=2.0) σ = sqrt(Ω / (K + 1) / 2) From 4ce14366d97f51bab90dea4cfddb7c5a85eb9cd6 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 19:48:52 +0800 Subject: [PATCH 03/24] Remove Rician inner constructor based on review suggestion Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 15b5affe9..69666aec5 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -25,7 +25,6 @@ External links: struct Rician{T<:Real} <: ContinuousUnivariateDistribution ν::T σ::T - Rician{T}(ν, σ) where {T} = new{T}(ν, σ) end function Rician(ν::T, σ::T; check_args=true) where {T<:Real} From e717bed3a1fa6d743fa405889f392dfb8388f478 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 19:53:04 +0800 Subject: [PATCH 04/24] Removed @inline as suggested during review Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 69666aec5..235e99aad 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -58,7 +58,7 @@ shape(d::Rician) = d.ν^2 / (2 * d.σ^2) scale(d::Rician) = d.ν^2 + 2 * d.σ^2 params(d::Rician) = (d.ν, d.σ) -@inline partype(d::Rician{T}) where {T<:Real} = T +partype(d::Rician{T}) where {T<:Real} = T #### Statistics From c0550871179c1b1a91025e71b3d05083cb47bfcd Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 19:54:45 +0800 Subject: [PATCH 05/24] Type improvement based on review Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 235e99aad..65ea7329b 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -69,7 +69,10 @@ mean(d::Rician{T}) where T = d.σ * √T(π/2) * L½(-d.ν^2/(2 * d.σ^2)) var(d::Rician) = 2 * d.σ^2 + d.ν^2 - (π * d.σ^2 / 2) * L½(-d.ν^2/(2 * d.σ^2))^2 quantile(d::Rician, q::Real) = quantile_newton(d, q, mean(d)) -mode(d::Rician{T}) where T = _minimize_gss(x -> -pdf(d, x), zero(T), 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) From cb43ef908768544bc6554897df6c93f696103dd3 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 19:55:26 +0800 Subject: [PATCH 06/24] Removed duplicate implementation for Rician pdf Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 65ea7329b..88bcccc04 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -101,8 +101,6 @@ function logpdf(d::Rician, x::Real) log.(x) - log(σ²) - (x.^2 + d.ν^2)/(2σ²) + log.(besseli.(zero(eltype(x)), x .* (d.ν / σ²))) end -pdf(d::Rician, x::Real) = exp.(logpdf(d, x)) - function cdf(d::Rician, x::Real) isnan(x) && return x x ≤ 0.0 && return zero(eltype(x)) From e7f6b63f45cf182a3e219dd1ce973d3770bc5857 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 19:57:01 +0800 Subject: [PATCH 07/24] Removed Rician keyword constructor --- src/univariate/continuous/rician.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 69666aec5..d15a36eea 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -7,10 +7,16 @@ The *Rician distribution* with parameters `ν` and `σ` has probability density 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)}}\\ +\nu = \\sigma\\sqrt{2K} +``` + ```julia Rician() # Rician distribution with parameters ν=0 and σ=1 Rician(ν, σ) # Rician distribution with parameters ν and σ -Rician(; K, Ω) # Rician distribution with shape K and scale Ω params(d) # Get the parameters, i.e. (ν, σ) shape(d) # Get the shape parameter K = ν²/2σ² @@ -34,12 +40,6 @@ end Rician(ν::Real, σ::Real) = Rician(promote(ν, σ)...) -function Rician(; K=0.0, Ω=2.0) - σ = sqrt(Ω / (K + 1) / 2) - ν = √2 * σ * sqrt(K) - Rician(ν, σ) -end - @distr_support Rician 0.0 Inf #### Conversions @@ -63,10 +63,10 @@ params(d::Rician) = (d.ν, d.σ) #### Statistics # helper -L½(x) = exp(x/2) * ((1-x) * besseli(zero(x), -x/2) - x * besseli(oneunit(x), -x/2)) +_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) * L½(-d.ν^2/(2 * d.σ^2)) -var(d::Rician) = 2 * d.σ^2 + d.ν^2 - (π * d.σ^2 / 2) * L½(-d.ν^2/(2 * d.σ^2))^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)) mode(d::Rician{T}) where T = _minimize_gss(x -> -pdf(d, x), zero(T), mean(d)) From 4344ea13cd21c2597ec520d1a7bed9a32044e380 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 20:49:56 +0800 Subject: [PATCH 08/24] Type stability improvments for Rician --- src/univariate/continuous/rician.jl | 32 ++++++++++++++--------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index df5ba3e9c..04d184e1a 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -10,8 +10,7 @@ f(x; \\nu, \\sigma) = \\frac{x}{\\sigma^2} \\exp\\left( \\frac{-(x^2 + \\nu^2)}{ If shape and scale parameters `K` and `Ω` are given instead, `ν` and `σ` may be computed from them: ```math -\\sigma = \\sqrt{\frac{\\Omega}{2(K + 1)}}\\ -\nu = \\sigma\\sqrt{2K} +\\sigma = \\sqrt{\\frac{\\Omega}{2(K + 1)}}, \\quad \\nu = \\sigma\\sqrt{2K} ``` ```julia @@ -31,6 +30,7 @@ External links: struct Rician{T<:Real} <: ContinuousUnivariateDistribution ν::T σ::T + Rician{T}(ν, σ) where {T} = new{T}(ν, σ) end function Rician(ν::T, σ::T; check_args=true) where {T<:Real} @@ -94,18 +94,18 @@ end #### PDF/CDF function logpdf(d::Rician, x::Real) - isnan(x) && return x - x ≤ 0.0 && return -convert(eltype(x), Inf) - isinf(x) && return -convert(eltype(x), Inf) + (x ≤ 0 || isinf(x)) && return oftype(float(x), -Inf) σ² = d.σ^2 - log.(x) - log(σ²) - (x.^2 + d.ν^2)/(2σ²) + log.(besseli.(zero(eltype(x)), x .* (d.ν / σ²))) + log(x) - log(σ²) - (x^2 + d.ν^2)/(2σ²) + log(besseli(zero(x), x .* (d.ν / σ²))) end function cdf(d::Rician, x::Real) - isnan(x) && return x - x ≤ 0.0 && return zero(eltype(x)) - isinf(x) && return one(eltype(x)) - first(quadgk(x -> pdf(d, x), zero(eltype(x)), x; rtol=1e-8)) + 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 @@ -124,20 +124,20 @@ function fit(::Type{<:Rician}, x::AbstractArray{T}; tol=1e-12, maxiters=500) whe μ₂ = var(x) r = μ₁ / √μ₂ if r < sqrt(π/(4-π)) - ν = zero(T) - σ = scale(fit(Rayleigh, x)) + ν = zero(float(T)) + σ = scale(fit(Rayleigh, float.(x))) else - ξ(θ) = 2 + θ^2 - T(π)/8 * exp(-θ^2 / 2) * ((2 + θ^2) * besseli(0, θ^2 / 4) + θ^2 * besseli(1, θ^2 / 4))^2 + ξ(θ) = 2 + θ^2 - π/8 * exp(-θ^2 / 2) * ((2 + θ^2) * besseli(0, θ^2 / 4) + θ^2 * besseli(1, θ^2 / 4))^2 g(θ) = sqrt(ξ(θ) * (1+r^2) - 2) - θ = one(T) + θ = 1.0 for j ∈ 1:maxiters θ⁻ = θ θ = g(θ) abs(θ - θ⁻) < tol && break end ξθ = ξ(θ) - σ = sqrt(μ₂ / ξθ) - ν = sqrt(μ₁^2 + (ξθ - 2) * σ^2) + σ = convert(float(T), sqrt(μ₂ / ξθ)) + ν = convert(float(T), sqrt(μ₁^2 + (ξθ - 2) * σ^2)) end Rician(ν, σ) end From ac63f7bc59bc37b211a622dbc9e7cfbb26ec0a09 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 21:00:41 +0800 Subject: [PATCH 09/24] Add back Rician Integer constructor to ensure tests pass --- src/univariate/continuous/rician.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 04d184e1a..bb47cb627 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -39,6 +39,7 @@ function Rician(ν::T, σ::T; check_args=true) where {T<:Real} end Rician(ν::Real, σ::Real) = Rician(promote(ν, σ)...) +Rician(ν::Integer, σ::Integer) = Rician(float(ν), float(σ)) @distr_support Rician 0.0 Inf From 31211639d5a4ff54320aedb72ec98f92d0805ef9 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 30 Aug 2021 21:07:16 +0800 Subject: [PATCH 10/24] Fixed Rician broken test cases --- src/univariate/continuous/rician.jl | 1 + test/rician.jl | 7 ------- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index bb47cb627..c343322e8 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -38,6 +38,7 @@ function Rician(ν::T, σ::T; check_args=true) where {T<:Real} return Rician{T}(ν, σ) end +Rician() = Rician(0.0, 1.0) Rician(ν::Real, σ::Real) = Rician(promote(ν, σ)...) Rician(ν::Integer, σ::Integer) = Rician(float(ν), float(σ)) diff --git a/test/rician.jl b/test/rician.jl index c75318432..6bd28e5d8 100644 --- a/test/rician.jl +++ b/test/rician.jl @@ -37,11 +37,4 @@ @test d1 isa Rician{Float64} @test params(d1) == (0.0, 1.0) - d1 = Rician(K=0.5, Ω=300.0) - @test d1 isa Rician{Float64} - @test params(d1)[1] ≈ 10.0 - @test params(d1)[2] ≈ 10.0 - @test shape(d1) ≈ 0.5 - @test scale(d1) ≈ 300.0 - end From c22230174f83e233923b9bccceadf20d808ac0f7 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Tue, 31 Aug 2021 23:54:39 +0800 Subject: [PATCH 11/24] Updated Rician implementation to use NoncentralChisq --- src/univariate/continuous/rician.jl | 44 ++++++++++++++++------ test/rician.jl | 57 +++++++++++++++++++++++++++++ test/truncate.jl | 3 +- 3 files changed, 92 insertions(+), 12 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index c343322e8..1a543584a 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -70,7 +70,6 @@ _Lhalf(x) = exp(x/2) * ((1-x) * besseli(zero(x), -x/2) - x * besseli(oneunit(x), 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) @@ -93,21 +92,44 @@ function _minimize_gss(f, a, b; tol=1e-12) (b + a) / 2 end -#### PDF/CDF +#### PDF/CDF/quantile delegated to NoncentralChisq + +# NOTE: +# 1.0ν and 1.0x used to support x::Float32, since the underlying NoncentralChisq delegates +# to StatsFuns(rmath.jl) which doesn't have support for Float32 + +function quantile(d::Rician, x::Real) + ν, σ = params(d) + return sqrt(quantile(NoncentralChisq(2, (1.0ν / σ)^2), 1.0x)) * σ +end + +function cquantile(d::Rician, x::Real) + ν, σ = params(d) + return sqrt(cquantile(NoncentralChisq(2, (1.0ν / σ)^2), 1.0x)) * σ +end + +function pdf(d::Rician, x::Real) + (x ≤ 0 || isinf(x)) && return zero(1.0x) + ν, σ = params(d) + return 2 * x / σ^2 * pdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) +end 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.ν / σ²))) + (x ≤ 0 || isinf(x)) && return oftype(1.0x, -Inf) + ν, σ = params(d) + return log(2 * x / σ^2) + logpdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) 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 + x ≤ 0 && return zero(1.0x) + ν, σ = params(d) + return cdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) +end + +function logcdf(d::Rician, x::Real) + x ≤ 0 && return oftype(1.0x, -Inf) + ν, σ = params(d) + return logcdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) end #### Sampling diff --git a/test/rician.jl b/test/rician.jl index 6bd28e5d8..1f3e5d1c1 100644 --- a/test/rician.jl +++ b/test/rician.jl @@ -18,6 +18,11 @@ @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) + d1 = Rician(10.0, 10.0) + @test median(d1) == quantile(d1, 0.5) + x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) + @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) + x = rand(Rician(5.0, 5.0), 100000) d1 = fit(Rician, x) @test d1 isa Rician{Float64} @@ -37,4 +42,56 @@ @test d1 isa Rician{Float64} @test params(d1) == (0.0, 1.0) + @test pdf(d1, -Inf) == 0.0 + @test pdf(d1, -1) == 0.0 + @test pdf(d1, Inf) == 0.0 + @test isnan(pdf(d1, NaN)) + + @test logpdf(d1, -Inf) == -Inf + @test logpdf(d1, -1) == -Inf + @test logpdf(d1, Inf) == -Inf + @test isnan(logpdf(d1, NaN)) + + @test cdf(d1, -Inf) == 0.0 + @test cdf(d1, -1) == 0.0 + @test cdf(d1, Inf) == 1.0 + @test isnan(cdf(d1, NaN)) + + @test logcdf(d1, -Inf) == -Inf + @test logcdf(d1, -1) == -Inf + @test logcdf(d1, Inf) == 0.0 + @test isnan(logcdf(d1, NaN)) + + @inferred pdf(d1, -Inf32) + @inferred pdf(d1, -1) + @inferred pdf(d1, 1.0) + @inferred pdf(d1, 1.0f0) + @inferred pdf(d1, 1) + @inferred pdf(d1, 1//2) + @inferred pdf(d1, Inf) + + @inferred logpdf(d1, -Inf32) + @inferred logpdf(d1, -1) + @inferred logpdf(d1, 1.0) + @inferred logpdf(d1, 1.0f0) + @inferred logpdf(d1, 1) + @inferred logpdf(d1, 1//2) + @inferred logpdf(d1, Inf) + + @inferred cdf(d1, -Inf32) + @inferred cdf(d1, -1) + @inferred cdf(d1, 1.0) + @inferred cdf(d1, 1.0f0) + @inferred cdf(d1, 1) + @inferred cdf(d1, 1//2) + @inferred cdf(d1, Inf) + + @inferred logcdf(d1, -Inf32) + @inferred logcdf(d1, -1) + @inferred logcdf(d1, 1.0) + @inferred logcdf(d1, 1.0f0) + @inferred logcdf(d1, 1) + @inferred logcdf(d1, 1//2) + @inferred logcdf(d1, Inf) + end diff --git a/test/truncate.jl b/test/truncate.jl index 12bfefc22..d521d8ab1 100644 --- a/test/truncate.jl +++ b/test/truncate.jl @@ -81,7 +81,8 @@ function verify_and_test(d::UnivariateDistribution, dct::Dict, n_tsamples::Int) if !(typeof(d) in [Distributions.Truncated{Distributions.NoncentralChisq{Float64},Distributions.Continuous, Float64}, Distributions.Truncated{Distributions.NoncentralF{Float64},Distributions.Continuous, Float64}, Distributions.Truncated{Distributions.NoncentralT{Float64},Distributions.Continuous, Float64}, - Distributions.Truncated{Distributions.StudentizedRange{Float64},Distributions.Continuous, Float64}]) + Distributions.Truncated{Distributions.StudentizedRange{Float64},Distributions.Continuous, Float64}, + Distributions.Truncated{Distributions.Rician{Float64},Distributions.Continuous, Float64}]) @test isapprox(logpdf(d, Dual(float(x))), lp, atol=sqrt(eps())) end # NOTE: this test is disabled as StatsFuns.jl doesn't have generic support for cdf() From 72bbe6d3bf0ea9f805f31f22f017d851ed8bc2da Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:25:40 +0800 Subject: [PATCH 12/24] Improve type stability in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 1a543584a..77a75f91e 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -109,9 +109,9 @@ function cquantile(d::Rician, x::Real) end function pdf(d::Rician, x::Real) - (x ≤ 0 || isinf(x)) && return zero(1.0x) ν, σ = params(d) - return 2 * x / σ^2 * pdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) + result = 2 * x / σ^2 * pdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + return x < 0 || isinf(x) ? zero(result) : result end function logpdf(d::Rician, x::Real) From 1fd4c6ecb0b789b51811d9ac066cefb2fcaa7625 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:25:57 +0800 Subject: [PATCH 13/24] Improve type stability in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 77a75f91e..55ae15785 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -115,9 +115,9 @@ function pdf(d::Rician, x::Real) end function logpdf(d::Rician, x::Real) - (x ≤ 0 || isinf(x)) && return oftype(1.0x, -Inf) ν, σ = params(d) - return log(2 * x / σ^2) + logpdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) + result = log(2 * x / σ^2) + logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + return x < 0 || isinf(x) ? oftype(result, -Inf) : result end function cdf(d::Rician, x::Real) From a17e2589636457466601a2894b8e01812585ac4c Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:26:32 +0800 Subject: [PATCH 14/24] Improve type stability in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 55ae15785..1be688cf2 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -121,9 +121,8 @@ function logpdf(d::Rician, x::Real) end function cdf(d::Rician, x::Real) - x ≤ 0 && return zero(1.0x) ν, σ = params(d) - return cdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) + return cdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) end function logcdf(d::Rician, x::Real) From 55b2c79e406d9c35f72645f379d6f37e3cb2bbd0 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:26:51 +0800 Subject: [PATCH 15/24] Improve type stability in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 1be688cf2..ce1eb1c77 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -126,9 +126,8 @@ function cdf(d::Rician, x::Real) end function logcdf(d::Rician, x::Real) - x ≤ 0 && return oftype(1.0x, -Inf) ν, σ = params(d) - return logcdf(NoncentralChisq(2, (1.0ν / σ)^2), (1.0x / σ)^2) + return logcdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) end #### Sampling From c1dced8f773717f0f36280dc7be25e0e77227f32 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:28:05 +0800 Subject: [PATCH 16/24] Avoid type conversion in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index ce1eb1c77..467118e79 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -147,7 +147,7 @@ function fit(::Type{<:Rician}, x::AbstractArray{T}; tol=1e-12, maxiters=500) whe r = μ₁ / √μ₂ if r < sqrt(π/(4-π)) ν = zero(float(T)) - σ = scale(fit(Rayleigh, float.(x))) + σ = scale(fit(Rayleigh, 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) From 397645dcb830c59dcbb50895069b39541231581a Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:29:25 +0800 Subject: [PATCH 17/24] Improve type stability in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 467118e79..b1c03ce92 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -151,8 +151,8 @@ function fit(::Type{<:Rician}, x::AbstractArray{T}; tol=1e-12, maxiters=500) whe 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(1) + for j in 1:maxiters θ⁻ = θ θ = g(θ) abs(θ - θ⁻) < tol && break From f8a7847586d533db33f587df53e4bb33a47bcf53 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Thu, 2 Sep 2021 01:50:16 +0800 Subject: [PATCH 18/24] Fixes edge cases without type instability in Rician --- src/univariate/continuous/rician.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index b1c03ce92..afab7c03f 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -94,10 +94,6 @@ end #### PDF/CDF/quantile delegated to NoncentralChisq -# NOTE: -# 1.0ν and 1.0x used to support x::Float32, since the underlying NoncentralChisq delegates -# to StatsFuns(rmath.jl) which doesn't have support for Float32 - function quantile(d::Rician, x::Real) ν, σ = params(d) return sqrt(quantile(NoncentralChisq(2, (1.0ν / σ)^2), 1.0x)) * σ @@ -116,18 +112,20 @@ end function logpdf(d::Rician, x::Real) ν, σ = params(d) - result = log(2 * x / σ^2) + logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + result = log(2 * abs(x) / σ^2) + logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) return x < 0 || isinf(x) ? oftype(result, -Inf) : result end function cdf(d::Rician, x::Real) ν, σ = params(d) - return cdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + result = cdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + return x < 0 ? zero(result) : result end function logcdf(d::Rician, x::Real) ν, σ = params(d) - return logcdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + result = logcdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2) + return x < 0 ? oftype(result, -Inf) : result end #### Sampling From 2ded077be7ba9c7b0fa099b35e27d1ea79dcbd12 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sun, 5 Sep 2021 14:04:38 +0800 Subject: [PATCH 19/24] =?UTF-8?q?Use=20StatsFuns.sqrthalf=CF=80=20in=20Ric?= =?UTF-8?q?ian?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index afab7c03f..a37f98cdc 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -67,7 +67,7 @@ partype(d::Rician{T}) where {T<:Real} = T # 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)) +mean(d::Rician) = d.σ * sqrthalfπ * _Lhalf(-d.ν^2/(2 * d.σ^2)) var(d::Rician) = 2 * d.σ^2 + d.ν^2 - (π * d.σ^2 / 2) * _Lhalf(-d.ν^2/(2 * d.σ^2))^2 function mode(d::Rician) From ebadc70af38053242bd24303c346cca6b74eee31 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sun, 5 Sep 2021 14:04:53 +0800 Subject: [PATCH 20/24] =?UTF-8?q?Use=20StatsFuns.half=CF=80=20in=20Rician?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index a37f98cdc..571362d93 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -68,7 +68,7 @@ partype(d::Rician{T}) where {T<:Real} = T _Lhalf(x) = exp(x/2) * ((1-x) * besseli(zero(x), -x/2) - x * besseli(oneunit(x), -x/2)) mean(d::Rician) = d.σ * sqrthalfπ * _Lhalf(-d.ν^2/(2 * d.σ^2)) -var(d::Rician) = 2 * d.σ^2 + d.ν^2 - (π * d.σ^2 / 2) * _Lhalf(-d.ν^2/(2 * d.σ^2))^2 +var(d::Rician) = 2 * d.σ^2 + d.ν^2 - halfπ * d.σ^2 * _Lhalf(-d.ν^2/(2 * d.σ^2))^2 function mode(d::Rician) m = mean(d) From 8bc514a3122678d71253590189810796626d33d6 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sun, 5 Sep 2021 14:05:27 +0800 Subject: [PATCH 21/24] Remove forcing of Float64 in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index 571362d93..c447ffc13 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -96,7 +96,7 @@ end function quantile(d::Rician, x::Real) ν, σ = params(d) - return sqrt(quantile(NoncentralChisq(2, (1.0ν / σ)^2), 1.0x)) * σ + return sqrt(quantile(NoncentralChisq(2, (ν / σ)^2), x)) * σ end function cquantile(d::Rician, x::Real) From dec17caf39421c503411021ed50a59d33877c7fe Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sun, 5 Sep 2021 14:05:44 +0800 Subject: [PATCH 22/24] Remove forcing of Float64 in Rician Co-authored-by: David Widmann --- src/univariate/continuous/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl index c447ffc13..9a5ccdb59 100644 --- a/src/univariate/continuous/rician.jl +++ b/src/univariate/continuous/rician.jl @@ -101,7 +101,7 @@ end function cquantile(d::Rician, x::Real) ν, σ = params(d) - return sqrt(cquantile(NoncentralChisq(2, (1.0ν / σ)^2), 1.0x)) * σ + return sqrt(cquantile(NoncentralChisq(2, (ν / σ)^2), x)) * σ end function pdf(d::Rician, x::Real) From 05a9699f4080b70b35f8d9e4caac1c494a815a71 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sun, 5 Sep 2021 14:14:13 +0800 Subject: [PATCH 23/24] Force same random numbers in Rician tests --- test/rician.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/rician.jl b/test/rician.jl index 1f3e5d1c1..a9e82913d 100644 --- a/test/rician.jl +++ b/test/rician.jl @@ -23,7 +23,7 @@ x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) - x = rand(Rician(5.0, 5.0), 100000) + x = rand(MersenneTwister(3456789), Rician(5.0, 5.0), 100000) d1 = fit(Rician, x) @test d1 isa Rician{Float64} @test params(d1)[1] ≈ 5.0 atol=1e-1 From 8a0cadacbc9c2a287b9dca812521c1994d12b168 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sun, 5 Sep 2021 16:13:49 +0800 Subject: [PATCH 24/24] Increase tolerence on Rician random tests --- test/rician.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/rician.jl b/test/rician.jl index a9e82913d..5674373b3 100644 --- a/test/rician.jl +++ b/test/rician.jl @@ -23,11 +23,11 @@ x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90]) - x = rand(MersenneTwister(3456789), Rician(5.0, 5.0), 100000) + 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 + @test params(d1)[1] ≈ 5.0 atol=0.2 + @test params(d1)[2] ≈ 5.0 atol=0.2 d1 = Rician(10.0f0, 10.0f0) @test d1 isa Rician{Float32}