From 0f02f93aa8f9dd419bbe2735284040b92ed236d4 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 14 Jul 2026 15:29:24 -0500 Subject: [PATCH 1/2] Add reliability diagnostics to GVar MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A GVar's moments are exact only when each operation is at most quadratic over the spread of its input and the input really is Gaussian. Neither holds in general. GVar now carries two fields that accumulate through arithmetic and say how far the result has drifted from a Gaussian: κ3 the third cumulant, exposed standardized as `skewness` err an estimate of the absolute error in `center`, exposed as `moment_error`, tracking the leading Taylor term each operation neglects `distrust` combines them into one dimensionless number: the displacement of the quantiles in units of `rad`. It is zero for an exact Gaussian and grows without bound as the description degrades. The two components are complementary. `exp` of a Gaussian is exactly lognormal, so its `moment_error` is zero while its skewness is large; a quadratic has neither. Two fixes fall out of the same derivation: Var[f(x)] = f'^2 σ^2 + f''^2 σ^4/2 + f' f''' σ^4 + O(σ^6). The third term was omitted although it is the same order as the second, and it is positive for `inv`, `log`, `sqrt` and `x^p`, so σ was systematically understated: `log(3 ± 0.5)` was low by 3.1%, now 0.5%. ThickNumbers builds the result of `hull`, `intersect`, `typemin` and `typemax` with the two-argument call `TN(lo, hi)`, but `GVar(center, σ)` is a different parametrization; those four now have GVar methods. The generic versions produced a hull spanning [-2.5, 6.5] for operands spanning [2, 4] and [2.5, 4.5], and an intersection wider than either operand. Integer powers use exact Gaussian raw moments, which terminate for a polynomial, rather than an asymptotic correction to the variance. Mean, variance and third cumulant are then all exact, and SpecialFunctions is no longer needed. Exponents above 20 overflow the moment sum and throw. Assisted-by: Claude Opus 4.8 --- Project.toml | 2 - README.md | 45 ++++- src/GaussianRandomVariables.jl | 340 ++++++++++++++++++++++++--------- test/runtests.jl | 140 +++++++++++--- 4 files changed, 401 insertions(+), 126 deletions(-) diff --git a/Project.toml b/Project.toml index c87206e..e86b328 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,10 @@ authors = ["Tim Holy and contributors"] version = "1.0.0-DEV" [deps] -SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" ThickNumbers = "b57aa878-5b76-4266-befc-f8e007760995" [compat] HypothesisTests = "0.11" -SpecialFunctions = "2" StableRNGs = "1" Statistics = "1" Test = "1" diff --git a/README.md b/README.md index 9fa4516..ce5442e 100644 --- a/README.md +++ b/README.md @@ -11,17 +11,54 @@ julia> x = 5 ± 1 5.0 ± 1.0 julia> 1/x -0.20800000000000002 ± 0.041569219381653054 +0.20800000000000002 ± 0.04595650117230423 (distrust 0.16) ``` -Related package: [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl) is a much more fully-developed and featureful package which also handles arithmetic with Gaussian random variables. However, it implements first-order (linear) error propagation, which leads to different rules of arithmetic: compare +## Reliability + +Arithmetic on `GVar`s is exact whenever each operation is at most quadratic over the spread of its input, and whenever the input really is Gaussian. Neither holds in general, and when they fail the reported `center` and `σ` can be badly wrong, most dangerously by making the value look more tightly determined than it is. + +A `GVar` therefore carries two diagnostics that accumulate through arithmetic: + +- `skewness(a)`, the standardized third cumulant, is zero for a true Gaussian. It measures how far the result has drifted from the shape a `GVar` can represent. Its sign is informative: a large negative skewness is a long tail toward *low* values. +- `moment_error(a)` estimates the absolute error in `center`, by tracking the leading Taylor term each operation neglects. + +`distrust(a)` combines them into one dimensionless number, roughly "how far `a`'s quantiles are displaced, in units of `rad(a)`". Smaller is better; zero means exactly Gaussian. + +```julia +julia> distrust((3 ± 0.5)^2 + 1) # quadratic: handled exactly +0.0 + +julia> distrust(exp(1 ± 0.1)) # gently curved over a narrow spread +0.05029... + +julia> distrust(exp(1 ± 1.5)) # wildly non-Gaussian; do not believe it +5.578... +``` + +This makes `GVar` usable for global optimization, where the mean and radius over a box are used to decide whether the box can contain the optimum. The moments alone will happily miss a narrow, deep minimum — a "slot canyon" — and prune a box that in fact holds the answer. `distrust` gives you a reason to subdivide anyway: + +```julia +y = f(lohi(GVar, lo, hi)) +if distrust(y) < 1 && mid(y) - 2*rad(y) > best_so_far + # safe to prune +else + # bisect regardless of what the parameters say +end +``` + +`distrust` is a heuristic, not a bound — it cannot make `GVar` rigorous. In particular it cannot see the error caused by treating repeated appearances of a variable as independent, as in `x*x` (write `x^2`). For guaranteed enclosures, use interval arithmetic. + +## Relation to other packages + +[Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl) is a much more fully-developed and featureful package which also handles arithmetic with Gaussian random variables. However, it implements first-order (linear) error propagation, which leads to different rules of arithmetic: compare ```julia julia> using GaussianRandomVariables julia> (0 ± 1)^2 -1.0 ± 1.4142135623730951 +1.0 ± 1.4142135623730951 (distrust 0.47) ``` with @@ -33,6 +70,8 @@ julia> (0 ± 1)^2 0.0 ± 0.0 ``` +The `GVar` answer is exact: the square of a standard normal is χ²₁, with mean 1 and standard deviation √2. Its distrust is nonzero because χ²₁ is strongly skewed (its skewness is √8), so `mid ∓ rad` is a poor description of where its values actually lie, and the distrust value alerts you to this. + Measurements is recommended for most users, but GaussianRandomVariables can be recommended if second-order accuracy matters in your application. GaussianRandomVariables is built on top of [ThickNumbers](https://github.com/timholy/ThickNumbers.jl), and the API for working with `GVar`s is described in detail there. diff --git a/src/GaussianRandomVariables.jl b/src/GaussianRandomVariables.jl index 93f354b..bbb0ade 100644 --- a/src/GaussianRandomVariables.jl +++ b/src/GaussianRandomVariables.jl @@ -1,39 +1,66 @@ module GaussianRandomVariables using ThickNumbers -using SpecialFunctions: SpecialFunctions, gamma import Base: +, -, *, /, //, ^, inv -import Base: abs, abs2, max, min, sin, cos, sincos, sqrt +import Base: abs, abs2, max, min, sqrt import Base: log, exp export GVar, ± +export skewness, moment_error, distrust const BaseReals = Union{AbstractFloat, Integer, AbstractIrrational, Rational} ## Type and constructors +""" + GVar(center, σ) + center ± σ + +A Gaussian random variable with mean `center` and standard deviation `σ`. + +Arithmetic on `GVar`s propagates the distribution of the result to second order, +so `(0 ± 1)^2` is `1 ± √2` rather than the first-order answer `0 ± 0`. + +A `GVar` additionally carries two *reliability diagnostics*, which accumulate +through arithmetic and say how far the result has drifted from being a genuine +Gaussian: + +- [`skewness`](@ref), the standardized third cumulant. Zero for a Gaussian. +- [`moment_error`](@ref), an estimate of the absolute error in `center`. + +[`distrust`](@ref) combines them into a single dimensionless number. When it is +large, `center` and `σ` do not describe the true distribution and should not be +believed. +""" struct GVar{T<:Real} <: ThickNumber{T} center::T σ::T + κ3::T # third cumulant, equal to the third central moment + err::T # estimated error in `center`, accumulated from Taylor truncation - function GVar{T}(center::Real, σ::Real) where T<:Real - # σ < zero(σ) && throw(DomainError(σ, "σ must be nonnegative")) - new{T}(center, σ) + function GVar{T}(center, σ, κ3, err) where T<:Real + new{T}(center, σ, κ3, err) end end -GVar{T}(x::Real) where T = GVar{T}(x, zero(x)) -GVar(x::T) where T<:Real = GVar{T}(x) +GVar{T}(center, σ) where T<:Real = GVar{T}(center, σ, zero(T), zero(T)) +GVar{T}(x::Real) where T<:Real = GVar{T}(x, zero(T)) +GVar(center::T, σ::T, κ3::T, err::T) where T<:Real = GVar{T}(center, σ, κ3, err) +GVar(center::Real, σ::Real, κ3::Real, err::Real) = GVar(promote(center, σ, κ3, err)...) GVar(center::T, σ::T) where T<:Real = GVar{T}(center, σ) GVar(center::Real, σ::Real) = GVar(promote(center, σ)...) +GVar(x::T) where T<:Real = GVar{T}(x) +# Integer and irrational parameters float-promote: the moment formulas divide. GVar(center::T, σ::T) where T<:Integer = GVar(float(center), float(σ)) GVar(center::T, σ::T) where T<:Irrational = GVar(float(center), float(σ)) +GVar(c::T, s::T, k::T, e::T) where T<:Integer = GVar(float(c), float(s), float(k), float(e)) +GVar(c::T, s::T, k::T, e::T) where T<:Irrational = GVar(float(c), float(s), float(k), float(e)) GVar(x::GVar) = x -GVar{T}(x::GVar) where T = GVar{T}(x.center, x.σ) +GVar{T}(x::GVar) where T = GVar{T}(x.center, x.σ, x.κ3, x.err) ±(center::Real, σ::Real) = GVar(center, σ) @@ -41,7 +68,64 @@ Base.promote_rule(::Type{GVar{T}}, ::Type{GVar{S}}) where {T<:Real,S<:Real} = GV Base.promote_rule(::Type{GVar{T}}, ::Type{S}) where {T<:Real,S<:BaseReals} = GVar{promote_type(T,S)} Base.AbstractFloat(a::GVar{<:AbstractFloat}) = a -Base.AbstractFloat(a::GVar) = GVar(AbstractFloat(a.center), AbstractFloat(a.σ)) +Base.AbstractFloat(a::GVar) = GVar(AbstractFloat(a.center), AbstractFloat(a.σ), + AbstractFloat(a.κ3), AbstractFloat(a.err)) + +## Reliability diagnostics + +""" + skewness(a::GVar) + +The standardized third cumulant of `a`, `κ3/σ^3`. It is zero for a true Gaussian +and grows as arithmetic drives the distribution away from one. + +A nonzero skewness means the span `mid(a) ± k*rad(a)` is misplaced: to leading +order the `k`-sigma quantile sits at `mid(a) + rad(a)*(k + (k^2-1)*skewness(a)/6)`. +A large negative skewness is a long tail toward *low* values, i.e. the objective +may dip far below `loval(a)`. +""" +function skewness(a::GVar{T}) where T + iszero(a.κ3) && return zero(float(T)) + return a.κ3 / a.σ^3 +end + +""" + moment_error(a::GVar) + +An estimate of the absolute error in `mid(a)`, accumulated across the operations +that produced `a`. + +Arithmetic on `GVar`s is exact when each operation is at most quadratic over the +spread of its input, so this tracks the leading *neglected* Taylor term, +`|f''''(c)|σ^4/8`, propagated forward by `|f'(c)|` at each subsequent step. +""" +moment_error(a::GVar) = a.err + +""" + distrust(a::GVar) + +A dimensionless measure of how much the Gaussian parameters of `a` can be +believed. It is zero when `a` is exactly Gaussian and grows without bound as the +description degrades; **smaller is better**. + +Roughly, `distrust` is the displacement of `a`'s quantiles measured in units of +`rad(a)`: values ≲0.01 mean `mid(a)` and `rad(a)` are solid, while values ≳1 mean +they are fiction and the domain must be subdivided regardless of what they say. + +It sums the two defects, which are also available on their own: +`moment_error(a)/rad(a)` (Taylor truncation) and `abs(skewness(a))/6` (the +leading Cornish–Fisher displacement from non-Gaussianity). + +`distrust` is a heuristic, not a bound. It cannot see the error caused by +treating repeated appearances of a variable as independent, as in `x*x`. +""" +function distrust(a::GVar{T}) where T + iszero(a.err) && iszero(a.κ3) && return zero(float(T)) + d = a.err / a.σ + abs(a.κ3) / (6 * a.σ^3) + # A NaN here means the moments carry no information at all (e.g. `σ` is + # infinite); that is the least trustworthy state, not an unremarkable one. + return isnan(d) ? convert(float(T), Inf) : d +end ## ThickNumbers API @@ -53,17 +137,37 @@ ThickNumbers.basetype(::Type{GVar{T}}) where T = GVar ThickNumbers.basetype(::Type{GVar}) = GVar ThickNumbers.emptyset(::Type{G}) where G<:GVar = G(0, -1) -## Display +# ThickNumbers' generic `hull`, `intersect`, `typemin` and `typemax` construct a +# result with the two-argument call `TN(lo, hi)`, but `GVar`'s two-argument +# constructor is `GVar(center, σ)`. Build these from the span directly. +Base.typemin(::Type{GVar{T}}) where T<:Real = GVar{T}(typemin(T), zero(T)) +Base.typemax(::Type{GVar{T}}) where T<:Real = GVar{T}(typemax(T), zero(T)) + +# A hull or intersection is a set operation, not a Gaussian: it has no meaningful +# skew, and it must not gain reliability that neither operand had. +_span(::Type{G}, lo, hi, err) where G<:GVar = G((lo + hi)/2, (hi - lo)/2, zero(lo), err) -Base.show(io::IO, a::GVar) = print(io, a.center, " ± ", a.σ) +function Base.intersect(a::GVar{T}, b::GVar{T}) where T + isdisjoint(a, b) && return emptyset(GVar{T}) + return _span(GVar{T}, max(loval(a), loval(b)), min(hival(a), hival(b)), max(a.err, b.err)) +end +ThickNumbers.hull(a::GVar{T}, b::GVar{T}) where T = + _span(GVar{T}, min(loval(a), loval(b)), max(hival(a), hival(b)), max(a.err, b.err)) +## Display + +function Base.show(io::IO, a::GVar) + print(io, a.center, " ± ", a.σ) + d = distrust(a) + iszero(d) || print(io, " (distrust ", round(d; sigdigits=2), ")") +end ## Trait functions and constants Base.zero(::GVar{T}) where T<:Real = zero(GVar{T}) Base.zero(::Type{GVar{T}}) where T<:Real = GVar(zero(T), zero(T)) Base.oneunit(::GVar{T}) where T<:Real = oneunit(GVar{T}) -Base.oneunit(::Type{GVar{T}}) where T<:Real = GVar(oneunit(T),zero(T)) +Base.oneunit(::Type{GVar{T}}) where T<:Real = GVar(oneunit(T), zero(T)) Base.real(a::GVar) = a Base.conj(a::GVar) = a @@ -71,36 +175,71 @@ Base.conj(a::GVar) = a function Base.hash(x::GVar, h::UInt) magic = Int === Int64 ? 0x21f1b1afbb07c31a : 0x237f8305 h += magic - return hash(x.center, hash(x.σ, h)) + return hash(x.center, hash(x.σ, hash(x.κ3, hash(x.err, h)))) end +## Moment propagation + +# For x ~ N(c, σ²) with third cumulant κ3, and smooth f, ordering terms by h where +# σ ~ h and κ3 ~ h³: +# +# E[f] = f + f''σ²/2 + f'''κ3/6 + O(h⁴) +# Var[f] = f'²σ² + f'f''κ3 + f''²σ⁴/2 + f'f'''σ⁴ + O(h⁵) +# κ3[f] = f'³κ3 + 3f'²f''σ⁴ + f''³σ⁶ + O(h⁷) +# +# The leading neglected term of the mean, f''''σ⁴/8, is not added to the result; +# it is accumulated into `err` as the error estimate. Both σ⁴ terms of the +# variance are the same order, so keeping one without the other understates σ. +function gmap(f::F, df::D1, ddf::D2, dddf::D3, d4f::D4, a::GVar) where {F,D1,D2,D3,D4} + c, σ, κ3, err = a.center, a.σ, a.κ3, a.err + f1, f2, f3, f4 = df(c), ddf(c), dddf(c), d4f(c) + σ2 = σ^2 + return assemble(f(c) + f2*σ2/2 + f3*κ3/6, + f1^2*σ2 + f1*f2*κ3 + f2^2*σ2^2/2 + f1*f3*σ2^2, + f1^3*κ3 + 3*f1^2*f2*σ2^2 + f2^3*σ2^3, + abs(f1)*err + abs(f4)*σ2^2/8) +end -## Arithmetic +# A negative variance means the expansion has broken down completely, which +# happens only when the input is far from Gaussian. Surrender the moments rather +# than report a plausible-looking answer. +function assemble(m, v, κ3, err) + v < zero(v) && return GVar(m, zero(v), zero(κ3), oftype(err, Inf)) + return GVar(m, sqrt(v), κ3, err) +end -ispos(x::Real) = x > zero(x) -isneg(x::Real) = x < zero(x) -isnonneg(x::Real) = x >= zero(x) -isnonpos(x::Real) = x <= zero(x) +# E[(c + δ)^n] for δ ~ N(0, σ²). Odd powers of δ vanish and E[δ^k] = σ^k (k-1)!!, +# so the sum terminates: for polynomial f the Gaussian moments are exact. +function gaussrawmoment(n::Integer, c::T, σ::T) where T<:AbstractFloat + s = zero(T) + binom = one(T) # C(n, k) + dfact = one(T) # (k-1)!! + for k in 0:2:n + if k > 0 + binom *= T(n - k + 2) * T(n - k + 1) / (T(k) * T(k - 1)) + dfact *= T(k - 1) + end + s += binom * c^(n - k) * σ^k * dfact + end + return s +end ## Addition and subtraction -function +(a::GVar{<:Real}, b::Real) - GVar(a.center + b, a.σ) -end + ++(a::GVar{<:Real}, b::Real) = GVar(a.center + b, a.σ, a.κ3, a.err) +(a::GVar{<:Integer}, b::Integer) = float(a) + float(b) -+(b::Real, a::GVar{<:Real}) = a+b ++(b::Real, a::GVar{<:Real}) = a + b -function -(a::GVar{<:Real}, b::Real) - GVar(a.center - b, a.σ) -end +-(a::GVar{<:Real}) = GVar(-a.center, a.σ, -a.κ3, a.err) +-(a::GVar{<:Real}, b::Real) = GVar(a.center - b, a.σ, a.κ3, a.err) -(a::GVar{<:Integer}, b::Integer) = float(a) - float(b) -function -(b::Real, a::GVar{<:Real}) - GVar(b - a.center, a.σ) -end +-(b::Real, a::GVar{<:Real}) = GVar(b - a.center, a.σ, -a.κ3, a.err) -(b::Integer, a::GVar{<:Integer}) = float(b) - float(a) +# Cumulants of independent variables add. function +(a::GVar{<:Real}, b::GVar{<:Real}) @fastmath begin - ret = GVar(a.center + b.center, sqrt(a.σ^2 + b.σ^2)) + ret = GVar(a.center + b.center, sqrt(a.σ^2 + b.σ^2), a.κ3 + b.κ3, a.err + b.err) return (isempty(a) | isempty(b)) ? emptyset(ret) : ret end end @@ -108,122 +247,137 @@ end function -(a::GVar{<:Real}, b::GVar{<:Real}) @fastmath begin - ret = GVar(a.center - b.center, sqrt(a.σ^2 + b.σ^2)) + ret = GVar(a.center - b.center, sqrt(a.σ^2 + b.σ^2), a.κ3 - b.κ3, a.err + b.err) return (isempty(a) | isempty(b)) ? emptyset(ret) : ret end end -(a::GVar{<:Integer}, b::GVar{<:Integer}) = float(a) - float(b) ## Multiplication + # Restrict to BaseReals so we avoid ambiguities with ForwardDiff -function *(x::BaseReals, a::GVar{<:Real}) - @fastmath begin - c, σ = x*a.center, abs(x)*a.σ - return GVar(c, σ) - end -end +*(x::BaseReals, a::GVar{<:Real}) = GVar(x*a.center, abs(x)*a.σ, x^3*a.κ3, abs(x)*a.err) *(a::GVar{<:Real}, x::BaseReals) = x*a # Prevent overflow by promoting to float *(x::Integer, a::GVar{<:Integer}) = float(x)*float(a) *(a::GVar{<:Integer}, x::Integer) = float(x)*float(a) +# For independent a and b these moments are exact given the first three cumulants +# of each: no fourth moment appears. function *(a::GVar{<:Real}, b::GVar{<:Real}) @fastmath begin - ret = GVar(a.center*b.center, sqrt(a.center^2*b.σ^2 + b.center^2*a.σ^2 + a.σ^2*b.σ^2)) + ca, σa, κa, ea = a.center, a.σ, a.κ3, a.err + cb, σb, κb, eb = b.center, b.σ, b.κ3, b.err + v = ca^2*σb^2 + cb^2*σa^2 + σa^2*σb^2 + κ3 = ca^3*κb + cb^3*κa + κa*κb + 3*ca*σa^2*κb + 3*cb*σb^2*κa + 6*ca*cb*σa^2*σb^2 + ret = GVar(ca*cb, sqrt(v), κ3, abs(cb)*ea + abs(ca)*eb + ea*eb) return (isempty(a) | isempty(b)) ? emptyset(ret) : ret end end *(a::GVar{<:Integer}, b::GVar{<:Integer}) = float(a)*float(b) # prevent overflow ## Division -function /(a::GVar{<:Real}, x::Real) - @fastmath begin - c, σ = a.center/x, a.σ / abs(x) - return GVar(c, σ) - end -end + +/(a::GVar{<:Real}, x::Real) = GVar(a.center/x, a.σ/abs(x), a.κ3/x^3, a.err/abs(x)) function inv(a::GVar{T}) where T<:AbstractFloat # must be AbstractFloat so typemax gives Inf z = zero(1/oneunit(T)) isempty(a) && return emptyset(basetype(a){typeof(z)}) - zero(T) ∈ a && return GVar(iszero(a.center) ? z : 1/a.center, typemax(T)) - # As of 2023-10-11, the explicit case for `var[X/Y]` on - # https://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables#Second_moment - # is badly wrong. Use the general case below. - return GVar(meanvar(inv, x -> -1/x^2, x -> 2/x^3, a.center, a.σ)...) + # Straddling zero, the reciprocal has no finite moments at all. + zero(T) ∈ a && return GVar(iszero(a.center) ? z : 1/a.center, typemax(T), zero(T), typemax(T)) + return gmap(inv, x -> -1/x^2, x -> 2/x^3, x -> -6/x^4, x -> 24/x^5, a) end inv(a::GVar{<:Real}) = inv(float(a)) - /(a::Real, b::GVar{<:Real}) = a*inv(b) - /(a::GVar{<:Real}, b::GVar{<:Real}) = a*inv(b) - //(a::GVar, b::GVar) = a / b # to deal with rationals ## Powers + Base.literal_pow(::typeof(^), x::GVar, ::Val{0}) = oneunit(x) Base.literal_pow(::typeof(^), x::GVar, ::Val{1}) = x -function Base.literal_pow(::typeof(^), x::GVar, ::Val{2}) - c2, σ2 = x.center^2, x.σ^2 - return GVar(c2 + σ2, sqrt(4*c2*σ2 + 2*σ2^2)*sign(x.σ)) -end Base.literal_pow(::typeof(^), x::GVar, ::Val{p}) where p = x^p -function ^(a::GVar, p::Integer) + +function ^(a::GVar{T}, p::Integer) where T<:AbstractFloat isempty(a) && return a + p == 0 && return oneunit(a) + p == 1 && return a p < 0 && return inv(a^(-p)) - cp2 = a.center ^ (p-2) - c2, σ2 = a.center^2, a.σ^2 - # If |c| ≫ σ, then these are accurate: - c′ = c2*cp2 + p*(p-1)*σ2*cp2/2 - σ′² = p^2*σ2*cp2^2*c2 - # But if |c| ≪ σ and p is even, the dominant term is from σ^p. So let's just add that in: - if iseven(p) - σ′² += (gamma(p + 0.5) * 2^p / sqrt(pi) - 1) * a.σ^(2p) # the -1 is from subtracting ^2 - end - return GVar(c′, sqrt(σ′²)) + 3p <= 60 || throw(DomainError(p, "exponents above 20 overflow the moment sum; rescale or use a float exponent")) + c, σ, κ3, err = a.center, a.σ, a.κ3, a.err + # x^p is a polynomial, so the Gaussian moments terminate and are exact. The + # only error is the input's own, plus the leading correction for its skew. + m = gaussrawmoment(p, c, σ) + m2 = gaussrawmoment(2p, c, σ) + m3 = gaussrawmoment(3p, c, σ) + f1 = p*c^(p-1) + f2 = p*(p-1)*c^(p-2) + f3 = p >= 3 ? p*(p-1)*(p-2)*c^(p-3) : zero(T) + return assemble(m + f3*κ3/6, + (m2 - m^2) + f1*f2*κ3, + (m3 - 3*m*m2 + 2*m^3) + f1^3*κ3, + abs(f1)*err) end +^(a::GVar{<:Real}, p::Integer) = float(a)^p function ^(a::GVar, p::Real) isinteger(p) && return a^Int(p) isempty(a) && return a - return GVar(meanvar(x->x^p, x->p*x^(p-1), x -> p*(p-1)*x^(p-2), a.center, a.σ)...) + return gmap(x -> x^p, + x -> p*x^(p-1), + x -> p*(p-1)*x^(p-2), + x -> p*(p-1)*(p-2)*x^(p-3), + x -> p*(p-1)*(p-2)*(p-3)*x^(p-4), a) end ## Functions -function meanvar(f::F, df::DF, ddf::DDF, c, σ) where {F, DF, DDF} - # As of 2023-10-11, the Wikipedia page - # https://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables#Second_moment - # is quite confused about this: it presents three different formulas. Fortunately, one of them is correct! - σ2 = σ^2 - c′ = f(c) + ddf(c) * σ2 / 2 - σ′² = df(c)^2 * σ2 + ddf(c)^2 * σ2^2 / 2 - return c′, sqrt(σ′²) -end - -abs(a::GVar) = GVar(abs(a.center), a.σ) -abs2(a::GVar) = a^2 - -min(a::GVar, b::GVar) = lohi(GVar, min(loval(a), loval(b)), min(hival(a), hival(b))) -max(a::GVar, b::GVar) = lohi(GVar, max(loval(a), loval(b)), max(hival(a), hival(b))) -min(a::Real, b::GVar) = lohi(GVar, min(a, loval(b)), min(a, hival(b))) -min(a::GVar, b::Real) = min(b, a) -max(a::Real, b::GVar) = lohi(GVar, max(a, loval(b)), max(a, hival(b))) -max(a::GVar, b::Real) = max(b, a) -function Base.exp(a::GVar{<:AbstractFloat}) +function abs(a::GVar) isempty(a) && return a - return GVar(exp(a.center + a.σ^2/2), exp(a.center + a.σ^2)*sqrt(1 - exp(-a.σ^2))) + # Away from zero `abs` merely reflects. Straddling zero the result is a folded + # normal, which is not Gaussian at all; the center is then wrong by O(σ). + straddles = zero(a.center) ∈ a + err = straddles ? a.err + a.σ : a.err + return GVar(abs(a.center), a.σ, sign(a.center)*a.κ3, err) end +abs2(a::GVar) = a^2 -function Base.log(a::GVar{<:AbstractFloat}) - isempty(a) && return a - return GVar(meanvar(log, inv, x->-1/x^2, a.center, a.σ)...) +# `min` and `max` of Gaussians are not Gaussian. Keep the span, report no skew, +# and inherit the larger center error. +for (f, lo, hi) in ((:min, :min, :min), (:max, :max, :max)) + @eval begin + function $f(a::GVar, b::GVar) + r = lohi(GVar, $lo(loval(a), loval(b)), $hi(hival(a), hival(b))) + return GVar(r.center, r.σ, zero(r.κ3), max(a.err, b.err)) + end + function $f(a::Real, b::GVar) + r = lohi(GVar, $lo(a, loval(b)), $hi(a, hival(b))) + return GVar(r.center, r.σ, zero(r.κ3), b.err) + end + $f(a::GVar, b::Real) = $f(b, a) + end end -function Base.sqrt(a::GVar{<:AbstractFloat}) +function exp(a::GVar{<:AbstractFloat}) isempty(a) && return a - return GVar(meanvar(sqrt, x -> 1/(2*sqrt(x)), x -> -1/(4 * sqrt(x^3)), a.center, a.σ)...) + c, σ, κ3, err = a.center, a.σ, a.κ3, a.err + σ2 = σ^2 + # The exponential of a Gaussian is exactly lognormal, so there is no + # truncation in σ: only the input's skew and center error propagate. + e2 = exp(σ2) + m = exp(c + σ2/2) + v = (e2 - 1) * exp(2c + σ2) + k = (e2 - 1)^2 * (e2 + 2) * exp(3c + 3σ2/2) + f1 = exp(c) # f' = f'' = f''' for exp + return assemble(m + f1*κ3/6, v + f1^2*κ3, k + f1^3*κ3, m*err) end +log(a::GVar{<:AbstractFloat}) = + isempty(a) ? a : gmap(log, x -> 1/x, x -> -1/x^2, x -> 2/x^3, x -> -6/x^4, a) + +sqrt(a::GVar{<:AbstractFloat}) = + isempty(a) ? a : gmap(sqrt, x -> 1/(2*sqrt(x)), x -> -1/(4*sqrt(x^3)), + x -> 3/(8*sqrt(x^5)), x -> -15/(16*sqrt(x^7)), a) + end # module diff --git a/test/runtests.jl b/test/runtests.jl index 1bf0b7f..4ce02f1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,36 +26,120 @@ function testscalar(f, μ, σ; n=1000, pval = 1e-5, filter=Returns(true), rtol=n return isapprox(mean(y), mid(fg); rtol) && isapprox(std(y), rad(fg); rtol) end +# Third central moment of a large sample, for checking `κ3` propagation. +function testskew(f, μ, σ; n=10^7, rtol=0.05, filter=Returns(true)) + x = Base.filter(filter, μ .+ σ .* randn(rng, n)) + y = f.(x) + m = mean(y) + return isapprox(mean((y .- m).^3), f(GVar(μ, σ)).κ3; rtol) +end + ispositive(x) = x > 0 @testset "GaussianRandomVariables.jl" begin - x = 3 ± 1 - e = GVar(0, -1) # empty - @test isempty(e) - @test isempty(e^1) - @test isempty(e^2) - @test isempty(^(e, 2)) # non-literal evaluation - @test x^2 ⩪ ^(x, 2) # non-literal - @test mid(x*x) < mid(x^2) - @test rad(x*x) < rad(x^2) - @test rad(GVar(0, 1)^2) ≈ sqrt(2) - @test mid(x - x) == 0 - @test rad(x - x) ≈ sqrt(2) - @test mid(1/x) ≈ 1/mid(x) + rad(x)^2/mid(x)^3 - @test rad(1/x)^2 ≈ rad(x)^2 / mid(x)^4 + 2 * rad(x)^4 / mid(x)^6 - @test sqrt(x) ⩪ exp(0.5 * log(x)) rtol=1e-2 - @test sqrt(x) ⩪ x^0.5 - - for μ in (2, 3, 10), σ in (0.1, 0.25, 1) - # `GVar` propagates moments to second order, which holds only while the spread - # is small compared to the distance from `μ` to a singularity or a sign change. - # At μ = 3σ the sampled distribution of `log(x)` already has a tail reaching - # tens of radii below `mid`, and these tests rightly reject it. - μ >= 4σ || continue - @test testscalar(x -> x^2, μ, σ) - @test testscalar(x -> x^1.8, μ, σ; filter=ispositive) - @test testscalar(exp, μ, σ; rtol=0.03, n=10^6) - @test testscalar(log, μ, σ; filter=ispositive) - @test testscalar(sqrt, μ, σ; filter=ispositive) + @testset "arithmetic" begin + x = 3 ± 1 + e = GVar(0, -1) # empty + @test isempty(e) + @test isempty(e^1) + @test isempty(e^2) + @test isempty(^(e, 2)) # non-literal evaluation + @test x^2 ⩪ ^(x, 2) # non-literal + @test mid(x*x) < mid(x^2) + @test rad(x*x) < rad(x^2) + @test rad(GVar(0, 1)^2) ≈ sqrt(2) + @test mid(x - x) == 0 + @test rad(x - x) ≈ sqrt(2) + @test mid(1/x) ≈ 1/mid(x) + rad(x)^2/mid(x)^3 + # Var[1/x] = σ²/c⁴ + 2σ⁴/c⁶ + f'f'''σ⁴, and f'f''' = 6/c⁶ for f = inv + @test rad(1/x)^2 ≈ rad(x)^2/mid(x)^4 + 8*rad(x)^4/mid(x)^6 + @test sqrt(x) ⩪ exp(0.5 * log(x)) rtol=1e-2 + @test sqrt(x) ⩪ x^0.5 + @test mid(-x) == -mid(x) && rad(-x) == rad(x) + @test x^0 ⩪ oneunit(x) + @test x^1 ⩪ x + @test 1/(1/x) ⩪ x rtol=0.1 + end + + @testset "distributions" begin + for μ in (2, 3, 10), σ in (0.1, 0.25, 1) + # `GVar` propagates moments to second order, which holds only while the spread + # is small compared to the distance from `μ` to a singularity or a sign change. + # At μ = 3σ the sampled distribution of `log(x)` already has a tail reaching + # tens of radii below `mid`, and these tests rightly reject it. + μ >= 4σ || continue + @test testscalar(x -> x^2, μ, σ) + @test testscalar(x -> x^1.8, μ, σ; filter=ispositive) + @test testscalar(exp, μ, σ; rtol=0.03, n=10^6) + @test testscalar(log, μ, σ; filter=ispositive) + @test testscalar(sqrt, μ, σ; filter=ispositive) + end + end + + # x^p is a polynomial, so the Gaussian moments terminate: mean, variance and + # third cumulant are all exact. + @testset "exact moments of integer powers" begin + c, σ = 3.0, 0.75 + y = GVar(c, σ)^2 + @test mid(y) ≈ c^2 + σ^2 + @test rad(y)^2 ≈ 4c^2*σ^2 + 2σ^4 + @test y.κ3 ≈ 24c^2*σ^4 + 8σ^6 + @test moment_error(y) == 0 # nothing is neglected + # x² of a zero-mean Gaussian is σ²·χ²₁, whose skewness is √8 + @test skewness(GVar(0, 1)^2) ≈ sqrt(8) + @test testskew(x -> x^3, 2.0, 0.5) + end + + @testset "reliability diagnostics" begin + # A box straight from a bisection is a plain Gaussian: nothing to distrust. + b = lohi(GVar, 0.0, 1.0) + @test distrust(b) == 0 + @test skewness(b) == 0 + @test moment_error(b) == 0 + + # Quadratics are handled exactly, so they must not be flagged... + @test moment_error((3 ± 0.5)^2 + 1) == 0 + # ...while a strongly-curved function over a wide box must be. + @test distrust(exp(1 ± 1.5)) > 1 + @test distrust(sqrt(3 ± 0.01)) < 0.01 + + # distrust decreases monotonically as the domain is subdivided, which is + # what makes "bisect until trustworthy" terminate. + ds = [distrust(1 / ((lohi(GVar, 0.5 - w, 0.5 + w) - 0.2)^2 + 0.05)) for w in (0.4, 0.2, 0.1, 0.05, 0.025)] + @test issorted(ds; rev=true) + @test last(ds) < 0.1 + + # exp of a Gaussian is exactly lognormal: the moments are right, but the + # result is badly skewed, so mid ∓ rad is not a sensible span. + y = exp(1 ± 0.7) + @test moment_error(y) == 0 + @test skewness(y) > 2 + @test distrust(y) > 0.4 + + # Straddling zero, 1/x has no finite moments at all. + @test distrust(1 / (0 ± 1)) == Inf + end + + # ThickNumbers' generic `hull`/`intersect`/`typemin`/`typemax` build a result + # with `TN(lo, hi)`, but `GVar(center, σ)` is a different parametrization. + @testset "set operations" begin + a, b = GVar(3.0, 1.0), GVar(3.5, 1.0) # spans [2,4] and [2.5,4.5] + @test loval(hull(a, b)) ≈ 2.0 && hival(hull(a, b)) ≈ 4.5 + @test loval(intersect(a, b)) ≈ 2.5 && hival(intersect(a, b)) ≈ 4.0 + @test isempty(intersect(GVar(0.0, 1.0), GVar(10.0, 1.0))) + @test loval(typemax(GVar{Float64})) == Inf + @test hival(typemin(GVar{Float64})) == -Inf + # Combining two values must never manufacture reliability neither had. + u = exp(1 ± 1.5) + @test moment_error(hull(u, GVar(1.0, 1.0))) >= moment_error(u) + end + + @testset "constructors" begin + @test GVar{Float64}(1, 0) === GVar(1.0, 0.0) + @test GVar(1, 2) === GVar(1.0, 2.0) + @test GVar(π, π) isa GVar{Float64} + @test GVar{Float32}(GVar(1.0, 2.0)) isa GVar{Float32} + @test valuetype(GVar(1.0, 2.0)) === Float64 + @test convert(GVar{Float64}, 3) ⩪ GVar(3.0, 0.0) end end From ce246b60b0cec21a59155ed358e126867acf1846 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 14 Jul 2026 16:05:03 -0500 Subject: [PATCH 2/2] Take typemin/typemax from ThickNumbers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ThickNumbers 1.1.2 builds generic results with `lohi` and `midrad` rather than with a two-argument `TN(lo, hi)` call, so `GVar`'s midpoint/σ parametrization no longer needs its own `typemin` and `typemax`. `hull` and `intersect` stay, for a different reason: a result built from the span alone carries the diagnostics of a fresh box, so the generic versions would report `distrust` of zero for the hull of two values that are not to be trusted at all. Assisted-by: Claude Opus 4.8 --- Project.toml | 2 +- src/GaussianRandomVariables.jl | 10 +++------- test/runtests.jl | 9 +++++---- 3 files changed, 9 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index e86b328..d45a368 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ HypothesisTests = "0.11" StableRNGs = "1" Statistics = "1" Test = "1" -ThickNumbers = "1" +ThickNumbers = "1.1.2" julia = "1.10" [extras] diff --git a/src/GaussianRandomVariables.jl b/src/GaussianRandomVariables.jl index bbb0ade..6ad5cdf 100644 --- a/src/GaussianRandomVariables.jl +++ b/src/GaussianRandomVariables.jl @@ -137,14 +137,10 @@ ThickNumbers.basetype(::Type{GVar{T}}) where T = GVar ThickNumbers.basetype(::Type{GVar}) = GVar ThickNumbers.emptyset(::Type{G}) where G<:GVar = G(0, -1) -# ThickNumbers' generic `hull`, `intersect`, `typemin` and `typemax` construct a -# result with the two-argument call `TN(lo, hi)`, but `GVar`'s two-argument -# constructor is `GVar(center, σ)`. Build these from the span directly. -Base.typemin(::Type{GVar{T}}) where T<:Real = GVar{T}(typemin(T), zero(T)) -Base.typemax(::Type{GVar{T}}) where T<:Real = GVar{T}(typemax(T), zero(T)) - # A hull or intersection is a set operation, not a Gaussian: it has no meaningful -# skew, and it must not gain reliability that neither operand had. +# skew. ThickNumbers builds its generic result from the span alone, which would +# reset the diagnostics to zero and so manufacture reliability that neither +# operand had; carry the larger center error across instead. _span(::Type{G}, lo, hi, err) where G<:GVar = G((lo + hi)/2, (hi - lo)/2, zero(lo), err) function Base.intersect(a::GVar{T}, b::GVar{T}) where T diff --git a/test/runtests.jl b/test/runtests.jl index 4ce02f1..d921fb8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -120,8 +120,6 @@ ispositive(x) = x > 0 @test distrust(1 / (0 ± 1)) == Inf end - # ThickNumbers' generic `hull`/`intersect`/`typemin`/`typemax` build a result - # with `TN(lo, hi)`, but `GVar(center, σ)` is a different parametrization. @testset "set operations" begin a, b = GVar(3.0, 1.0), GVar(3.5, 1.0) # spans [2,4] and [2.5,4.5] @test loval(hull(a, b)) ≈ 2.0 && hival(hull(a, b)) ≈ 4.5 @@ -129,9 +127,12 @@ ispositive(x) = x > 0 @test isempty(intersect(GVar(0.0, 1.0), GVar(10.0, 1.0))) @test loval(typemax(GVar{Float64})) == Inf @test hival(typemin(GVar{Float64})) == -Inf - # Combining two values must never manufacture reliability neither had. - u = exp(1 ± 1.5) + # A result built from the span alone would report the diagnostics of a + # fresh box; combining values must not manufacture reliability. + u = 1 / ((3 ± 0.9)^2 - 8) # moment_error > 0 + @test moment_error(u) > 0 @test moment_error(hull(u, GVar(1.0, 1.0))) >= moment_error(u) + @test moment_error(intersect(u, GVar(mid(u), 10rad(u)))) >= moment_error(u) end @testset "constructors" begin