FastRicianLikelihoods

Rice(ν, σ)

The Rice distribution with shape parameter $\nu$ and scale parameter $\sigma$.

The probability density function is

\[p(x \mid \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)\]

where $I_0$ is the modified Bessel function of the first kind of order zero, and $x \ge 0$.

External links

• Rice distribution Wiki

besseli0(x::T) where {T <: Union{Float32, Float64}}

Modified Bessel function of the first kind of order zero, $I_0(x)$.

besseli0x(x::T) where {T <: Union{Float32, Float64}}

Scaled modified Bessel function of the first kind of order zero, $I_0(x) e^{-|x|}$.

besseli1(x::T) where {T <: Union{Float32, Float64}}

Modified Bessel function of the first kind of order one, $I_1(x)$.

besseli1i0(x::T) where {T <: Union{Float32, Float64}}

Ratio of modified Bessel functions of the first kind of orders one and zero, $I_1(x) / I_0(x)$.

besseli1i0m1(x::T) where {T <: Union{Float32, Float64}}

Ratio of modified Bessel functions minus one, $I_1(x) / I_0(x) - 1$.

besseli1i0x(x::T) where {T <: Union{Float32, Float64}}

Ratio of modified Bessel functions divided by argument, $(I_1(x) / I_0(x)) / x$.

besseli1x(x::T) where {T <: Union{Float32, Float64}}

Scaled modified Bessel function of the first kind of order one, $I_1(x) e^{-|x|}$.

besseli2(x::T) where {T <: Union{Float32, Float64}}

Modified Bessel function of the first kind of order two, $I_2(x)$.

besseli2x(x::T) where {T <: Union{Float32, Float64}}

Scaled modified Bessel function of the first kind of order two, $I_2(x) e^{-|x|}$.

logbesseli0(x::T) where {T <: Union{Float32, Float64}}

Log of modified Bessel function of the first kind of order zero, $\log I_0(x)$.

logbesseli0x(x::T) where T <: Union{Float32, Float64}

Log of scaled modified Bessel function of the first kind of order zero, $\log(I_0(x) e^{-|x|})$.

logbesseli1(x::T) where {T <: Union{Float32, Float64}}

Log of modified Bessel function of the first kind of order one, $\log I_1(x)$.

logbesseli1x(x::T) where {T <: Union{Float32, Float64}}

Log of scaled modified Bessel function of the first kind of order one, $\log(I_1(x) e^{-|x|})$.

logbesseli2(x::T) where {T <: Union{Float32, Float64}}

Log of modified Bessel function of the first kind of order two, $\log I_2(x)$.

logbesseli2x(x::T) where {T <: Union{Float32, Float64}}

Log of scaled modified Bessel function of the first kind of order two, $\log(I_2(x) e^{-|x|})$.

neglogpdf_qrician(x::Real, ν::Real, logσ::Real, δ::Real, order::Val)
neglogpdf_qrician(n::Int, ν::Real, logσ::Real, δ::Real, order::Val)
neglogpdf_qrician(x::Real, ν::Real, δ::Real, order::Val)

Negative log-probability mass function of the quantized Rician distribution.

Five-parameter form (real-valued argument)

For $\sigma = \exp(\log\sigma)$ and bin width $\delta$, the pmf is

\[p_{\mathrm{QRice}}(x \mid \nu, \sigma, \delta) = \int_{x}^{x+\delta} p_{\mathrm{Rice}}(y \mid \nu, \sigma) \, dy\]

Computes $-\log p_{\mathrm{QRice}}(x \mid \nu, \sigma, \delta)$ using $N$-point Gauss–Legendre quadrature with order::Val{N} where $N \ge 1$; the case $N = 1$ reduces to the midpoint rule.

Four-parameter form (unit scale)

The four-argument method sets $\sigma = 1$ and computes $-\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$.

Five-parameter form (discrete argument)

For integer argument n::Int, computes the negative log-probability at $x = n\delta$; equivalent to neglogpdf_qrician(n * δ, ν, logσ, δ, order).

See neglogpdf_rician.

neglogpdf_rician(x::Real, ν::Real, logσ::Real)
neglogpdf_rician(x::Real, ν::Real)

Negative log-density of the Rician distribution.

Three-parameter form

For $\sigma = \exp(\log\sigma)$, computes the negative log-density $-\log p_{\mathrm{Rice}}(x \mid \nu, \sigma)$:

\[p_{\mathrm{Rice}}(x \mid \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)\]

where $I_0$ is the modified Bessel function of the first kind of order zero, and $x \ge 0$.

Two-parameter form (unit scale)

The two-argument method sets $\sigma = 1$ and computes $f(x, \nu) = -\log p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1)$:

\[f(x,\nu) \coloneqq -\log p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1) = \frac{x^2+\nu^2}{2} - \log x - \log I_0(x\nu)\]

pdf_rician(x::Real, ν::Real, logσ::Real)
pdf_rician(x::Real, ν::Real)

Probability density function of the Rician distribution.

Three-parameter form

For $\sigma = \exp(\log\sigma)$, computes $p_{\mathrm{Rice}}(x \mid \nu, \sigma)$ as defined in neglogpdf_rician(x, ν, logσ).

Two-parameter form (unit scale)

The two-argument method sets $\sigma = 1$ and computes $p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1) = \exp(-f(x, \nu))$.

See neglogpdf_rician.

∇neglogpdf_qrician(x::Real, ν::Real, δ::Real, order::Val)

Gradient of the unit-scale quantized negative log-probability.

Computes $g = \nabla \Omega = (\Omega_x, \Omega_\nu, \Omega_\delta)$ where $\Omega = -\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$.

See neglogpdf_qrician.

∇neglogpdf_qrician_with_primal(x::Real, ν::Real, δ::Real, order::Val)

Primal value and gradient of the unit-scale quantized negative log-probability.

Computes $(\Omega, g)$ where

• $\Omega = -\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$,
• $g = \nabla \Omega = (\Omega_x, \Omega_\nu, \Omega_\delta)$.

See neglogpdf_qrician.

∇neglogpdf_rician(x::Real, ν::Real)

Gradient of the unit-scale negative log-density.

Computes $g = (f_x, f_\nu)$ where $f = -\log p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1)$.

See neglogpdf_rician.

∇pdf_rician(x::Real, ν::Real)

Gradient of the unit-scale Rician density with respect to $(x, \nu)$.

Computes $(p_x, p_\nu)$ where $p = p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1)$.

See pdf_rician, neglogpdf_rician.

∇²neglogpdf_qrician(x::Real, ν::Real, δ::Real, order::Val)

Hessian of the unit-scale quantized negative log-probability.

Computes $H = \mathrm{vech}(\nabla^2 \Omega) = (\Omega_{xx}, \Omega_{x\nu}, \Omega_{x\delta}, \Omega_{\nu\nu}, \Omega_{\nu\delta}, \Omega_{\delta\delta})$ where $\Omega = -\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$.

See neglogpdf_qrician.

∇²neglogpdf_qrician_with_gradient(x::Real, ν::Real, δ::Real, order::Val)

Gradient and Hessian of the unit-scale quantized negative log-probability.

Computes $(g, H)$ where

• $g = \nabla \Omega = (\Omega_x, \Omega_\nu, \Omega_\delta)$,
• $H = \mathrm{vech}(\nabla^2 \Omega) = (\Omega_{xx}, \Omega_{x\nu}, \Omega_{x\delta}, \Omega_{\nu\nu}, \Omega_{\nu\delta}, \Omega_{\delta\delta})$.

See neglogpdf_qrician.

∇²neglogpdf_qrician_with_primal_and_gradient(x::Real, ν::Real, δ::Real, order::Val)

Primal value, gradient, and Hessian of the unit-scale quantized negative log-probability.

Computes $(\Omega, g, H)$ where

• $\Omega = -\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$,
• $g = \nabla \Omega = (\Omega_x, \Omega_\nu, \Omega_\delta)$,
• $H = \mathrm{vech}(\nabla^2 \Omega) = (\Omega_{xx}, \Omega_{x\nu}, \Omega_{x\delta}, \Omega_{\nu\nu}, \Omega_{\nu\delta}, \Omega_{\delta\delta})$.

See neglogpdf_qrician.

∇²neglogpdf_rician(x::Real, ν::Real)

Hessian of the unit-scale negative log-density.

Computes $H = (f_{xx}, f_{x\nu}, f_{\nu\nu})$ where $f = -\log p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1)$.

See neglogpdf_rician.

∇²neglogpdf_rician_with_gradient(x::Real, ν::Real)

Gradient and Hessian of the unit-scale negative log-density.

Computes $(g, H)$ where

• $f = -\log p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1)$
• $g = (f_x, f_\nu)$,
• $H = (f_{xx}, f_{x\nu}, f_{\nu\nu})$,

See neglogpdf_rician.

∇³neglogpdf_qrician_jacobian_with_primal_gradient_and_hessian(x::Real, ν::Real, δ::Real, order::Val)

Jacobian of third-order derivatives of the unit-scale quantized negative log-probability.

Computes $(\Omega, g, H, J)$ where

• $\Omega = -\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$,
• $g = \nabla \Omega = (\Omega_x, \Omega_\nu, \Omega_\delta)$,
• $H = \mathrm{vech}(\nabla^2 \Omega) = (\Omega_{xx}, \Omega_{x\nu}, \Omega_{x\delta}, \Omega_{\nu\nu}, \Omega_{\nu\delta}, \Omega_{\delta\delta})$,
• $J = \nabla H = \begin{bmatrix} \Omega_{xxx} & \Omega_{xx\nu} & \Omega_{xx\delta} \\ \Omega_{x\nu x} & \Omega_{x\nu\nu} & \Omega_{x\nu\delta} \\ \Omega_{x\delta x} & \Omega_{x\delta\nu} & \Omega_{x\delta\delta} \\ \Omega_{\nu\nu x} & \Omega_{\nu\nu\nu} & \Omega_{\nu\nu\delta} \\ \Omega_{\nu\delta x} & \Omega_{\nu\delta\nu} & \Omega_{\nu\delta\delta} \\ \Omega_{\delta\delta x} & \Omega_{\delta\delta\nu} & \Omega_{\delta\delta\delta} \end{bmatrix} \in \mathbb{R}^{6\times 3}$.

See neglogpdf_qrician.

∇³neglogpdf_qrician_vjp_with_primal_gradient_and_hessian(Δ::SVector{6, <:Real}, x::Real, ν::Real, δ::Real, order::Val)

Vector-Jacobian product for third-order derivatives of the unit-scale quantized negative log-probability.

Computes $(\Omega, g, H, J^T \Delta)$ where

• $\Omega = -\log p_{\mathrm{QRice}}(x \mid \nu, \sigma = 1, \delta)$,
• $g = \nabla \Omega = (\Omega_x, \Omega_\nu, \Omega_\delta)$,
• $H = \mathrm{vech}(\nabla^2 \Omega) = (\Omega_{xx}, \Omega_{x\nu}, \Omega_{x\delta}, \Omega_{\nu\nu}, \Omega_{\nu\delta}, \Omega_{\delta\delta})$,
• $J^T \Delta = (\nabla H)^T \Delta = \begin{bmatrix} \Omega_{xxx} & \Omega_{x\nu x} & \Omega_{x\delta x} & \Omega_{\nu\nu x} & \Omega_{\nu\delta x} & \Omega_{\delta\delta x} \\ \Omega_{xx\nu} & \Omega_{x\nu\nu} & \Omega_{x\delta\nu} & \Omega_{\nu\nu\nu} & \Omega_{\nu\delta\nu} & \Omega_{\delta\delta\nu} \\ \Omega_{xx\delta} & \Omega_{x\nu\delta} & \Omega_{x\delta\delta} & \Omega_{\nu\nu\delta} & \Omega_{\nu\delta\delta} & \Omega_{\delta\delta\delta} \end{bmatrix} \Delta \in \mathbb{R}^3$.

See neglogpdf_qrician.

∇³neglogpdf_rician_with_gradient_and_hessian(x::Real, ν::Real)

Gradient, Hessian, and third-order partial derivatives of the unit-scale negative log-density.

Computes $(g, H, T)$ where

• $f = -\log p_{\mathrm{Rice}}(x \mid \nu, \sigma = 1)$.
• $g = (f_x, f_\nu)$,
• $H = (f_{xx}, f_{x\nu}, f_{\nu\nu})$,
• $T = (f_{xxx}, f_{xx\nu}, f_{x\nu\nu}, f_{\nu\nu\nu})$,

See neglogpdf_rician.

GaussHalfHermite

Compute nodes x and weights w for Gauss–Half–Hermite quadrature on [0, ∞) with weight w(x) = x^γ * exp(-x^2).

\[\int_{0}^{\infty} x^{\gamma} e^{-x^{2}} f(x) \, dx \approx \sum_{i=1}^{N} w_i f(x_i)\]

Numerical method:

• Stable recurrence coefficients (α_n, β_n) via Ball’s reparameterization g_n and Newton’s method (tridiagonal Jacobian; O(N) per iteration).
• Nodes/weights from the symmetric Jacobi matrix via Golub–Welsch.

Public API:

• gausshalfhermite_gw(N, γ; normalize=false) -> x, w
• gausshalfhermite_rec_coeffs(N, γ) -> α, β

References:

• Ball J. (2002) SIAM J. Numer. Anal. 40:2311–2317.
• Golub GH, Welsch JH. (1969) Math. Comp. 23:221–230.
• Shizgal B. (1981) J. Comput. Phys. 41:309–328.
• Galant D. (1969) Math. Comp. 23:674–s39.

gausshalfhermite_gw(N, γ; normalize = false) -> (x, w)

Nodes x and weights w for N‑point Gauss–Half–Hermite quadrature.

\[\int_{0}^{\infty} x^{\gamma} e^{-x^{2}} f(x) \, dx \approx \sum_{i=1}^{N} w_i f(x_i)\]

Uses the Golub–Welsch algorithm on the symmetric Jacobi matrix from (α, β) coefficients computed by gausshalfhermite_rec_coeffs.

If normalize=true, rescale weights and nodes to correspond to weighting function w(x) = x^γ * exp(-x^2 / 2) / √(2π) by setting x ← √2 * x and w ← 2^{γ/2} w / √π.

Arguments:

• N::Integer
• γ::Real (γ > -1)

Keyword arguments:

• normalize::Bool=false

Returns:

• (x, w): nodes and weights

Example

julia> x, w = gausshalfhermite_gw(2, 2.0);

julia> f(x) = x^2;

julia> I = dot(w, f.(x));

julia> I ≈ 3 * sqrt(π) / 8 # ∫_{0}^{∞} x^4 * exp(-x^2) dx
true

gausshalfhermite_rec_coeffs(N, γ) -> (α, β)

Recurrence coefficients for monic polynomials orthogonal w.r.t. w(x) = x^γ * exp(-x^2) on [0, ∞) using the three-term recurrence formula P_{n+1}(x) = (x - α_n) P_n(x) - β_n P_{n-1}(x).

Arguments:

• N::Integer: number of coefficients; returns α₀:α_{N-1} and β₀:β_{N-1}
• γ::Real: exponent in the weight (γ > -1)

Returns:

• (α, β): diagonal α and off-diagonal squares β of the Jacobi matrix

GaussLegendre

Compute nodes x and weights w for Gauss–Legendre quadrature on [-1, 1].

\[\int_{-1}^{1} f(x)\,dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Numerical method:

• For n ≤ 60, roots are found by Newton’s method; Legendre polynomials are evaluated via a three‑term recurrence.
• For n > 60, an O(n) asymptotic expansion provides nodes and weights; optional refinement by Newton is available via refine=true.

Public API:

• gausslegendre(n::Integer, ::Type{T} = Float64; refine = true) -> x, w

Porting notes:

• Adapted from FastGaussQuadrature.jl (gausslegendre.jl)
• Generalized from Float64 to arbitrary T.
• Added refine keyword to control Newton refinement for large n.

gausslegendre(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Legendre quadrature.

\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Example

julia> x, w = gausslegendre(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 2/5 # ∫_{-1}^{1} x^4 dx
true