Skip to content

Commit 597190b

Browse files
feat(uq): add entropy, mgf, and cf methods for LineParametersPDF
1 parent fc14583 commit 597190b

1 file changed

Lines changed: 125 additions & 0 deletions

File tree

src/uq/distributions.jl

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -374,3 +374,128 @@ function Distributions.modes(d::LineParametersPDF{T}) where {T}
374374

375375
return [(d.edges[i] + d.edges[i+1]) / 2 for i in indices]
376376
end
377+
378+
"""
379+
entropy(d::LineParametersPDF)
380+
381+
Calculate the differential entropy (base e).
382+
H(X) = - ∫ f(x) * log(f(x)) dx
383+
= - Σ ∫_{e_i}^{e_{i+1}} [d_i * log(d_i)] dx
384+
= - Σ [d_i * log(d_i) * w_i]
385+
= - Σ [p_i * log(d_i)]
386+
where p_i = d_i * w_i is the probability mass of bin i.
387+
"""
388+
function Distributions.entropy(d::LineParametersPDF{T}) where {T}
389+
acc = zero(T)
390+
widths = diff(d.edges)
391+
392+
@inbounds for i in 1:length(d.dens)
393+
dens_i = d.dens[i]
394+
395+
# If density is 0, (f(x) * log(f(x))) -> 0.
396+
# So we just skip the bin.
397+
if dens_i > 0
398+
width_i = widths[i]
399+
prob_mass_i = dens_i * width_i
400+
401+
# This is base e (natural log)
402+
acc -= prob_mass_i * log(dens_i)
403+
end
404+
end
405+
return acc
406+
end
407+
408+
"""
409+
entropy(d::LineParametersPDF, b::Real)
410+
411+
Calculate the differential entropy with a specified base b.
412+
H_b(X) = H_e(X) / log(b)
413+
"""
414+
function Distributions.entropy(d::LineParametersPDF, b::Real)
415+
(b > 0 && b != 1) ||
416+
throw(ArgumentError("Entropy base must satisfy b > 0 and b ≠ 1, got b = $b"))
417+
418+
# Just convert the base e entropy.
419+
return Distributions.entropy(d) / log(b)
420+
end
421+
422+
423+
# ─────────────────────────────────────────────────────────────────────────────
424+
# Moment Generating & Characteristic Functions
425+
# ─────────────────────────────────────────────────────────────────────────────
426+
427+
"""
428+
mgf(d::LineParametersPDF, t::Real)
429+
430+
Moment Generating Function
431+
M(t) = E[e^(tX)] = ∫ e^(tx) * f(x) dx
432+
= Σ ∫_{e_i}^{e_{i+1}} [d_i * e^(tx)] dx
433+
= Σ d_i * [e^(tx) / t]_{e_i}^{e_{i+1}}
434+
= Σ d_i/t * (e^(t*e_{i+1}) - e^(t*e_i))
435+
436+
This is numerically catastrophic for t -> 0.
437+
We rewrite:
438+
Term_i = d_i * e^(t*e_i) * (e^(t*(e_{i+1}-e_i)) - 1) / t
439+
Let w_i = e_{i+1} - e_i.
440+
Term_i = d_i * e^(t*e_i) * w_i * (e^(t*w_i) - 1) / (t*w_i)
441+
Term_i = d_i * e^(t*e_i) * w_i * exprel(t*w_i)
442+
443+
`Base.Math.exprel(x)` is the numerically stable (e^x - 1) / x.
444+
"""
445+
function Distributions.mgf(d::LineParametersPDF{T}, t::Real) where {T}
446+
t == 0 && return one(T) # MGF(0) = 1
447+
448+
acc = zero(T)
449+
450+
@inbounds for i in 1:length(d.dens)
451+
dens_i = d.dens[i]
452+
dens_i == 0 && continue
453+
454+
a = d.edges[i]
455+
b = d.edges[i+1]
456+
w = b - a
457+
458+
# This is the argument to exprel: z = t*w
459+
z = t * w
460+
461+
# Calculate (e^z - 1) / z, handling z=0
462+
# This is the stable implementation
463+
exprel_z = iszero(z) ? one(z) : expm1(z) / z
464+
465+
# Term_i = d_i * e^(t*a) * w_i * exprel(t*w_i)
466+
acc += dens_i * exp(t * a) * w * exprel_z
467+
end
468+
return acc
469+
end
470+
471+
472+
"""
473+
cf(d::LineParametersPDF, t::Real)
474+
475+
Characteristic Function
476+
ϕ(t) = E[e^(itX)] = MGF(it)
477+
478+
This is the exact same derivation as the MGF, just
479+
substituting `t` with `it`.
480+
"""
481+
function Distributions.cf(d::LineParametersPDF{T}, t::Real) where {T}
482+
t == 0 && return one(Complex{T})
483+
484+
acc = zero(Complex{T})
485+
486+
@inbounds for i in 1:length(d.dens)
487+
dens_i = d.dens[i]
488+
dens_i == 0 && continue
489+
490+
a = d.edges[i]
491+
b = d.edges[i+1]
492+
w = b - a
493+
494+
z = im * t * w
495+
496+
exprel_z = iszero(z) ? one(z) : expm1(z) / z
497+
498+
acc += dens_i * exp(im * t * a) * w * exprel_z
499+
end
500+
return acc
501+
end

0 commit comments

Comments
 (0)