|
1 | 1 | module GaussianRandomVariables |
2 | 2 |
|
3 | | -# Write your package code here. |
| 3 | +using ThickNumbers |
| 4 | +using SpecialFunctions: SpecialFunctions, gamma |
4 | 5 |
|
| 6 | +import Base: +, -, *, /, //, ^, inv |
| 7 | +import Base: abs, abs2, max, min, sin, cos, sincos, sqrt |
| 8 | +import Base: log, exp |
| 9 | + |
| 10 | +export GVar, ± |
| 11 | + |
| 12 | +const BaseReals = Union{AbstractFloat, Integer, AbstractIrrational, Rational} |
| 13 | + |
| 14 | +## Type and constructors |
| 15 | + |
| 16 | +struct GVar{T<:Real} <: ThickNumber{T} |
| 17 | + center::T |
| 18 | + σ::T |
| 19 | + |
| 20 | + function GVar{T}(center::Real, σ::Real) where T<:Real |
| 21 | + # σ < zero(σ) && throw(DomainError(σ, "σ must be nonnegative")) |
| 22 | + new{T}(center, σ) |
| 23 | + end |
| 24 | +end |
| 25 | + |
| 26 | +GVar{T}(x::Real) where T = GVar{T}(x, zero(x)) |
| 27 | +GVar(x::T) where T<:Real = GVar{T}(x) |
| 28 | + |
| 29 | +GVar(center::T, σ::T) where T<:Real = GVar{T}(center, σ) |
| 30 | +GVar(center::Real, σ::Real) = GVar(promote(center, σ)...) |
| 31 | + |
| 32 | +GVar(center::T, σ::T) where T<:Integer = GVar(float(center), float(σ)) |
| 33 | +GVar(center::T, σ::T) where T<:Irrational = GVar(float(center), float(σ)) |
| 34 | + |
| 35 | +GVar(x::GVar) = x |
| 36 | +GVar{T}(x::GVar) where T = GVar{T}(x.center, x.σ) |
| 37 | + |
| 38 | +±(center::Real, σ::Real) = GVar(center, σ) |
| 39 | + |
| 40 | +Base.promote_rule(::Type{GVar{T}}, ::Type{GVar{S}}) where {T<:Real,S<:Real} = GVar{promote_type(T,S)} |
| 41 | +Base.promote_rule(::Type{GVar{T}}, ::Type{S}) where {T<:Real,S<:BaseReals} = GVar{promote_type(T,S)} |
| 42 | + |
| 43 | +AbstractFloat(a::GVar{<:AbstractFloat}) = a |
| 44 | +AbstractFloat(a::GVar) = GVar(AbstractFloat(a.center), AbstractFloat(a.σ)) |
| 45 | + |
| 46 | +## ThickNumbers API |
| 47 | + |
| 48 | +ThickNumbers.loval(a::GVar) = a.center - a.σ |
| 49 | +ThickNumbers.hival(a::GVar) = a.center + a.σ |
| 50 | +ThickNumbers.lohi(::Type{G}, lo, hi) where G<:GVar = G((lo + hi)/2, nextfloat((hi - lo)/2)) |
| 51 | +ThickNumbers.midrad(::Type{G}, center, σ) where G<:GVar = G(center, σ) |
| 52 | +ThickNumbers.basetype(::Type{GVar{T}}) where T = GVar |
| 53 | +ThickNumbers.basetype(::Type{GVar}) = GVar |
| 54 | +ThickNumbers.emptyset(::Type{G}) where G<:GVar = G(0, -1) |
| 55 | + |
| 56 | +## Display |
| 57 | + |
| 58 | +Base.show(io::IO, a::GVar) = print(io, a.center, " ± ", a.σ) |
| 59 | + |
| 60 | + |
| 61 | +## Trait functions and constants |
| 62 | + |
| 63 | +Base.zero(::GVar{T}) where T<:Real = zero(GVar{T}) |
| 64 | +Base.zero(::Type{GVar{T}}) where T<:Real = GVar(zero(T), zero(T)) |
| 65 | +Base.oneunit(::GVar{T}) where T<:Real = oneunit(GVar{T}) |
| 66 | +Base.oneunit(::Type{GVar{T}}) where T<:Real = GVar(oneunit(T),zero(T)) |
| 67 | + |
| 68 | +Base.real(a::GVar) = a |
| 69 | +Base.conj(a::GVar) = a |
| 70 | + |
| 71 | +function Base.hash(x::GVar, h::UInt) |
| 72 | + magic = Int === Int64 ? 0x21f1b1afbb07c31a : 0x237f8305 |
| 73 | + h += magic |
| 74 | + return hash(x.center, hash(x.σ, h)) |
| 75 | +end |
| 76 | + |
| 77 | + |
| 78 | +## Arithmetic |
| 79 | + |
| 80 | +ispos(x::Real) = x > zero(x) |
| 81 | +isneg(x::Real) = x < zero(x) |
| 82 | +isnonneg(x::Real) = x >= zero(x) |
| 83 | +isnonpos(x::Real) = x <= zero(x) |
| 84 | + |
| 85 | +## Addition and subtraction |
| 86 | +function +(a::GVar{<:Real}, b::Real) |
| 87 | + GVar(a.center + b, a.σ) |
| 88 | +end |
| 89 | ++(a::GVar{<:Integer}, b::Integer) = float(a) + float(b) |
| 90 | ++(b::Real, a::GVar{<:Real}) = a+b |
| 91 | + |
| 92 | +function -(a::GVar{<:Real}, b::Real) |
| 93 | + GVar(a.center - b, a.σ) |
| 94 | +end |
| 95 | +-(a::GVar{<:Integer}, b::Integer) = float(a) - float(b) |
| 96 | +function -(b::Real, a::GVar{<:Real}) |
| 97 | + GVar(b - a.center, a.σ) |
| 98 | +end |
| 99 | +-(b::Integer, a::GVar{<:Integer}) = float(b) - float(a) |
| 100 | + |
| 101 | +function +(a::GVar{<:Real}, b::GVar{<:Real}) |
| 102 | + @fastmath begin |
| 103 | + ret = GVar(a.center + b.center, sqrt(a.σ^2 + b.σ^2)) |
| 104 | + return (isempty(a) | isempty(b)) ? emptyset(ret) : ret |
| 105 | + end |
| 106 | +end |
| 107 | ++(a::GVar{<:Integer}, b::GVar{<:Integer}) = float(a) + float(b) |
| 108 | + |
| 109 | +function -(a::GVar{<:Real}, b::GVar{<:Real}) |
| 110 | + @fastmath begin |
| 111 | + ret = GVar(a.center - b.center, sqrt(a.σ^2 + b.σ^2)) |
| 112 | + return (isempty(a) | isempty(b)) ? emptyset(ret) : ret |
| 113 | + end |
| 114 | +end |
| 115 | +-(a::GVar{<:Integer}, b::GVar{<:Integer}) = float(a) - float(b) |
| 116 | + |
| 117 | +## Multiplication |
| 118 | +# Restrict to BaseReals so we avoid ambiguities with ForwardDiff |
| 119 | +function *(x::BaseReals, a::GVar{<:Real}) |
| 120 | + @fastmath begin |
| 121 | + c, σ = x*a.center, abs(x)*a.σ |
| 122 | + return GVar(c, σ) |
| 123 | + end |
| 124 | +end |
| 125 | +*(a::GVar{<:Real}, x::BaseReals) = x*a |
| 126 | +# Prevent overflow by promoting to float |
| 127 | +*(x::Integer, a::GVar{<:Integer}) = float(x)*float(a) |
| 128 | +*(a::GVar{<:Integer}, x::Integer) = float(x)*float(a) |
| 129 | + |
| 130 | +function *(a::GVar{<:Real}, b::GVar{<:Real}) |
| 131 | + @fastmath begin |
| 132 | + ret = GVar(a.center*b.center, sqrt(a.center^2*b.σ^2 + b.center^2*a.σ^2 + a.σ^2*b.σ^2)) |
| 133 | + return (isempty(a) | isempty(b)) ? emptyset(ret) : ret |
| 134 | + end |
| 135 | +end |
| 136 | +*(a::GVar{<:Integer}, b::GVar{<:Integer}) = float(a)*float(b) # prevent overflow |
| 137 | + |
| 138 | +## Division |
| 139 | +function /(a::GVar{<:Real}, x::Real) |
| 140 | + @fastmath begin |
| 141 | + c, σ = a.center/x, a.σ / abs(x) |
| 142 | + return GVar(c, σ) |
| 143 | + end |
| 144 | +end |
| 145 | + |
| 146 | +function inv(a::GVar{T}) where T<:AbstractFloat # must be AbstractFloat so typemax gives Inf |
| 147 | + z = zero(1/oneunit(T)) |
| 148 | + isempty(a) && return emptyset(basetype(a){typeof(z)}) |
| 149 | + zero(T) ∈ a && return GVar(iszero(a.center) ? z : 1/a.center, typemax(T)) |
| 150 | + # As of 2023-10-11, the explicit case for `var[X/Y]` on |
| 151 | + # https://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables#Second_moment |
| 152 | + # is badly wrong. Use the general case below. |
| 153 | + return GVar(meanvar(inv, x -> -1/x^2, x -> 2/x^3, a.center, a.σ)...) |
5 | 154 | end |
| 155 | +inv(a::GVar{<:Real}) = inv(float(a)) |
| 156 | + |
| 157 | + |
| 158 | +/(a::Real, b::GVar{<:Real}) = a*inv(b) |
| 159 | + |
| 160 | +/(a::GVar{<:Real}, b::GVar{<:Real}) = a*inv(b) |
| 161 | + |
| 162 | +//(a::GVar, b::GVar) = a / b # to deal with rationals |
| 163 | + |
| 164 | +## Powers |
| 165 | +Base.literal_pow(::typeof(^), x::GVar, ::Val{0}) = oneunit(x) |
| 166 | +Base.literal_pow(::typeof(^), x::GVar, ::Val{1}) = x |
| 167 | +function Base.literal_pow(::typeof(^), x::GVar, ::Val{2}) |
| 168 | + c2, σ2 = x.center^2, x.σ^2 |
| 169 | + return GVar(c2 + σ2, sqrt(4*c2*σ2 + 2*σ2^2)*sign(x.σ)) |
| 170 | +end |
| 171 | +Base.literal_pow(::typeof(^), x::GVar, ::Val{p}) where p = x^p |
| 172 | +function ^(a::GVar, p::Integer) |
| 173 | + isempty(a) && return a |
| 174 | + p < 0 && return inv(a^(-p)) |
| 175 | + cp2 = a.center ^ (p-2) |
| 176 | + c2, σ2 = a.center^2, a.σ^2 |
| 177 | + # If |c| ≫ σ, then these are accurate: |
| 178 | + c′ = c2*cp2 + p*(p-1)*σ2*cp2/2 |
| 179 | + σ′² = p^2*σ2*cp2^2*c2 |
| 180 | + # But if |c| ≪ σ and p is even, the dominant term is from σ^p. So let's just add that in: |
| 181 | + if iseven(p) |
| 182 | + σ′² += (gamma(p + 0.5) * 2^p / sqrt(pi) - 1) * a.σ^(2p) # the -1 is from subtracting <a^p>^2 |
| 183 | + end |
| 184 | + return GVar(c′, sqrt(σ′²)) |
| 185 | +end |
| 186 | + |
| 187 | +function ^(a::GVar, p::Real) |
| 188 | + isinteger(p) && return a^Int(p) |
| 189 | + isempty(a) && return a |
| 190 | + return GVar(meanvar(x->x^p, x->p*x^(p-1), x -> p*(p-1)*x^(p-2), a.center, a.σ)...) |
| 191 | +end |
| 192 | + |
| 193 | +## Functions |
| 194 | +function meanvar(f::F, df::DF, ddf::DDF, c, σ) where {F, DF, DDF} |
| 195 | + # As of 2023-10-11, the Wikipedia page |
| 196 | + # https://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables#Second_moment |
| 197 | + # is quite confused about this: it presents three different formulas. Fortunately, one of them is correct! |
| 198 | + σ2 = σ^2 |
| 199 | + c′ = f(c) + ddf(c) * σ2 / 2 |
| 200 | + σ′² = df(c)^2 * σ2 + ddf(c)^2 * σ2^2 / 2 |
| 201 | + return c′, sqrt(σ′²) |
| 202 | +end |
| 203 | + |
| 204 | +abs(a::GVar) = GVar(abs(a.center), a.σ) |
| 205 | +abs2(a::GVar) = a^2 |
| 206 | + |
| 207 | +min(a::GVar, b::GVar) = lohi(GVar, min(loval(a), loval(b)), min(hival(a), hival(b))) |
| 208 | +max(a::GVar, b::GVar) = lohi(GVar, max(loval(a), loval(b)), max(hival(a), hival(b))) |
| 209 | +min(a::Real, b::GVar) = lohi(GVar, min(a, loval(b)), min(a, hival(b))) |
| 210 | +min(a::GVar, b::Real) = min(b, a) |
| 211 | +max(a::Real, b::GVar) = lohi(GVar, max(a, loval(b)), max(a, hival(b))) |
| 212 | +max(a::GVar, b::Real) = max(b, a) |
| 213 | + |
| 214 | +function Base.exp(a::GVar{<:AbstractFloat}) |
| 215 | + isempty(a) && return a |
| 216 | + return GVar(meanvar(exp, exp, exp, a.center, a.σ)...) |
| 217 | +end |
| 218 | +Base.exp(a::GVar) = exp(float(a)) |
| 219 | + |
| 220 | +function Base.log(a::GVar{<:AbstractFloat}) |
| 221 | + isempty(a) && return a |
| 222 | + return GVar(meanvar(log, inv, x->-1/x^2, a.center, a.σ)...) |
| 223 | +end |
| 224 | + |
| 225 | +function Base.sqrt(a::GVar{<:AbstractFloat}) |
| 226 | + isempty(a) && return a |
| 227 | + return GVar(meanvar(sqrt, x -> 1/(2*sqrt(x)), x -> -1/(4 * sqrt(x^3)), a.center, a.σ)...) |
| 228 | +end |
| 229 | + |
| 230 | +end # module |
0 commit comments