Repository navigation
Evaluate fdistlogpdf symmetrically in the two degrees of freedom - #226
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #226 +/- ##
==========================================
+ Coverage 75.50% 75.91% +0.40%
==========================================
Files 23 23
Lines 1135 1146 +11
==========================================
+ Hits 857 870 +13
+ Misses 278 276 -2 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
The density was evaluated as
(xlogy(ν1, ν1 / ν2) + xlogy(ν1 - 2, x)
- xlogy(ν1 + ν2, 1 + (ν1 / ν2) * x)) / 2 - logbeta(ν1 / 2, ν2 / 2)
which breaks down at both ends:
* for large ν2 the term (ν1 / ν2) * x is tiny and the explicit 1 +
discards it - the argument rounds to exactly 1.0 from ν2 ≈ 1e17 on,
so the factor drops out. fdistlogpdf(1, 1e18, 6.6) returned -1.862
instead of -5.162, i.e. the density was off by e^3.3 ≈ 27.
* for large ν1 the first and last terms both grow like ν1 * log(ν1)
and cancel. fdistlogpdf(1e16, 5, 2) returned -5.914 instead of
-1.670.
Writing the density in terms of the beta variate u = ν1x / (ν1x + ν2),
with log(u) = -log1p(ν2 / (ν1x)) and log1p(-u) = -log1p(ν1x / ν2), is
symmetric in ν1 and ν2 and removes both. The relative error stays below
1e-14 in either limit.
As a side effect fdistlogpdf(3, 7, Inf) now returns -Inf rather than NaN.
The x < 0 result moved from a shared line onto its own branch, which
exposed the fact that it was never executed: before this commit there
were no fdist tests below the support, and the ternary counted as covered
whenever x > 0 ran. A second testitem pins the boundary values - the sign
of x, ±Inf, NaN, x == 0 and the oftype conversion that keeps the
out-of-support branch from widening to Float64.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
4dbacb7 to
f8bbaed
Compare
andreasnoack
left a comment
andreasnoack
left a comment
There was a problem hiding this comment.
Thanks for this! I checked out the branch and exercised the new fdistlogpdf: the algebra matches the textbook form (a·log u + b·log1p(-u) - log x - logbeta(a, b)), all the new assertions pass, type stability holds for Int, Rational{Int}, Float16, Float32, Float64 and BigFloat on all three branches, NaN propagates correctly, and the existing ChainRulesCore rule already used the same log1p form, so primal and derivative agree. A 20k-point random sweep (ν1, ν2 ∈ [1e-3, 1e9], x ∈ [1e-6, 1e6]) against a 512-bit reference showed no accuracy regression apart from two ~5e-13 cases in the both-dof-large regime that #225 already tracks.
Two comments:
1. ν2 / (ν1 * x) overflows for small x
In the general branch, ν2 / (ν1 * x) overflows once x < ν2 / (ν1 * floatmax(T)), xlog1py then receives Inf, and the whole log-density collapses to -Inf. The previous formulation had no such term and was correct in this regime, so this is a regression:
| call | this PR | correct |
|---|---|---|
fdistlogpdf(1.0, 1e18, 1e-291) |
-Inf |
333.0 |
fdistlogpdf(3.0, 7.0, 1e-310) |
-Inf |
-356.07 |
fdistlogpdf(1f0, 7f0, 1f-38) |
-Inf |
42.79 |
fdistlogpdf(Float16(1), Float16(7), Float16(1e-4)) |
-Inf |
3.65 |
The Float16 case is the one I'd worry about most: x = 1e-4 is an entirely ordinary normal-range argument, but 7 / 1e-4 already exceeds floatmax(Float16) = 65504, so pdf silently returns zero probability. Note also that this bites inside the large-ν2 regime the PR targets — at ν1 = 1, ν2 = 1e18, x = 1e-290 still gives 332.96, and then x = 1e-291 falls off a cliff to -Inf.
One way out is to derive both terms from the single ratio r = ν1 * x / ν2 that the second term already needs:
r = ν1 * x / ν2
aterm = r > 1 ? -xlog1py(a, inv(r)) : xlogy(a, r) - xlog1py(a, r)Both branches are cancellation-free, and inv(r) cannot overflow.
test/accuracy.jl only probes x ∈ {0.5, 2.0, 6.6}, so it can't currently see this — worth adding a few tiny-x cases.
2. rtol = 1.0e-13 in test/accuracy.jl is quite tight
The tightest of the new assertions, (ν1, ν2, x) = (1e12, 5.0, 0.5), comes in at 2.9e-14 relative (and (1e18, 5.0, 0.5) at 2.4e-14), so there's only ~3.5x headroom. The error is dominated by logbeta, so a future SpecialFunctions change or a different platform's log could flip this to a CI failure. 1e-12 would still fail every one of these cases on master (errors there range from 1e-8 to 6e1), so it loses no regression-detecting power.
Very minor: the PR description says the accuracy item is 30 assertions; it's 60.
— reviewed and drafted by Claude Code on my behalf
…D rule `ν2 / (ν1 * x)` overflows for `x < ν2 / (ν1 * floatmax(T))`, which made the log-density collapse to `-Inf` - for `Float16(1e-4)`, an entirely ordinary argument. `log(u)` is now taken as `-log1p(1 / r)` only for `r = ν1 * x / ν2 > 1`; below that the previous textbook form is kept, with `log1p(r)` in place of `log(1 + r)`. That form also stays accurate when `r` underflows, and it absorbs the `x == 0` and `NaN` cases, so one branch disappears. The same ratio is used to build the beta variate of the cdf family and in the `∂/∂ν1` of the ChainRules rule, with the same overflow: on master `fdistlogcdf(Float16(1), Float16(7), Float16(1e-4))` is `-Inf` instead of `-4.87`. Both now split on the same halves of the domain. The tolerance of the accuracy tests is relaxed from `1e-13` to `1e-12`, which still fails every one of those assertions on master (their errors range from 1e-8 to 6e1) but leaves more than the previous 3.5x of headroom over the tightest case, whose error is dominated by `logbeta`. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Thanks, both confirmed and fixed in 555ec9a. 1. Small 2. While checking where else the ratio is formed: Assertion count corrected in the description (it is 143 now). One thing I have deliberately left alone: — drafted by Claude Code on my behalf |
|
You are right about I also confirmed the rest of Two things the new commit regresses, though, both in the opposite direction from the one it repairs — huge arguments rather than tiny ones. Both are correct on master:
1.
|
Three problems, all of the same kind as the ones already addressed on this
branch, and all found while re-reviewing it.
* the ratio `ν1 * x / ν2` overflows once `ν1 * x` does, even when the ratio
itself is perfectly representable, and the log-density then collapses to
`-Inf`. That is a regression against master, which formed the ratio as
`(ν1 / ν2) * x`: `fdistlogpdf(Float16(4), Float16(100), Float16(6e4))` was
`-Inf` instead of `-392.2`, for an entirely ordinary `Float16` argument.
The ratio is built from `ν1 / ν2` again, which the textbook branch needs
anyway, and `log1p(r)` is taken as `log(ν1 / ν2) + log(x)` where `r` itself
still overflows - both summands are positive whenever it does. That also
repairs `fdistlogpdf(3.0, 7.0, 1e308)`, `-Inf` on master as well.
* the branch condition `r > 1` was tighter than necessary. The symmetric form
only requires that `1 / r` not overflow, and widening the condition to
`isfinite(inv(r))` hands the whole ridge `x ≈ ν2 / ν1` - where the textbook
form distributes `ν1 / 2` over two logarithms that cancel - to the
well-conditioned branch. The relative error there drops from 1.3e-14 to
3.0e-16. The textbook branch is still reached for `x == 0`, for an
underflowing `r`, and for `NaN`.
* the beta variate of the cdf family overflowed in `y + ν2` for
`ν2 > floatmax / 2`, returning `u = 0`: `fdistcdf(1.0, 1e308, 1e308)` was
`0.0` instead of `0.5`. It is now built from the same ratio as the density.
In the `ChainRulesCore` rule, `∂/∂x` was a difference of two terms both growing
like `ν1 / x`, which cancelled completely - at `(1e18, 5, 2)` it returned `0.0`
where the derivative is `-1.125`. Combining them over `ν1 * x + ν2` collapses
the two into `-(ν1 * ν2 * a / 2 + 1) / x`, since `a * ν2 - b * ν1 * x` is
`(ν1 * ν2 / 2) * (1 - x)`; `a` is already in the setup. The relative error at
that point goes from 1 to 1.4e-18.
The reason to leave this alone was that the rewrite changes the value at
`x = Inf`, but that is a property of `a = (x - 1) / temp1`, which is `NaN`
there and wrong once `temp1` overflows, not of the rewrite. With `a` guarded,
`x = Inf` returns the correct limits for all three partials rather than `NaN`
for two of them, and `∂/∂ν1` and `∂/∂ν2` become accurate for a large `x` as
well.
The tests move out of `test/accuracy.jl` and into the files that already cover
this ground: the `BigFloat` references join `Utils` in `test/utils.jl`, the
accuracy assertions join the other `numerical issue` testsets in
`test/misc.jl`, and the boundary values become a testset of `Misc` next to
`support of binomial distribution`. The grid gains the ridge and the large-`x`
end, the three branches are checked for type stability across mixed argument
types, and the assertions on `fdistpdf` are dropped, since it is defined as
`exp(fdistlogpdf)`.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
In the cdf an overflowing `ν1 * x` only saturates the beta variate to 1, which is the right answer there, whereas in `fdistlogpdf` it collapses the density to `-Inf`. Dividing afterwards therefore costs nothing, and it keeps `u` exact where `ν1 / ν2` would go subnormal: `fdistcdf(1.0, floatmax(), floatmax())` builds `u = 0.5` rather than one ulp below it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`-(ν1 * ν2 * a / 2 + 1) / x` multiplies the two degrees of freedom directly, which overflows for degrees of freedom of `floatmax` scale even where the derivative is ordinary: at `(1e200, 1e300, 1e200)` it returned `-Inf` for a true `-5e99`. Grouping the product as `ν1 * (ν2 * a)` keeps every intermediate in range, since `ν2 * a` is a ratio of the two rather than a product, and it is the subexpression `∂/∂ν1` already computes. The grid now reaches `floatmax`-scale degrees of freedom in both directions, which is what exposed this. That needed the `BigFloat` references to evaluate at 1024 bits rather than the default 256: they are the textbook cancelling forms on purpose, and `logbeta(a, b)` cancels `lgamma(b) - lgamma(a + b)` at the hundredth digit when the two degrees of freedom are 100 orders of magnitude apart, while `∂/∂x` cancels to 300 digits. At 256 bits the reference for `logpdf(1e200, 1e300, 0.5)` is off by a relative 3.6e3 and the one for `dlogpdf_dx(1e300, 1e200, 2.0)` by a relative 1.0, so either would have validated anything. The pdf grid is 106 assertions now, 63 of which fail on master. The pairs are `1e200` against `1e300` rather than `floatmax` against `1e300`, where the `Float64` `logbeta` alone costs enough to reach 4.9e-13 - within a factor of two of the tolerance. The two `fdistcdf` values reported in review are pinned as well. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
1. The product. Your diagnosis was right, the remedy differs. Rather than recovering 2. The cdf branch. Also fixed, though split on 3. The grid. Done, and it earned its keep twice over. It caught a third problem, in the
And it caught a trap in the tests themselves: the The pdf grid is 106 assertions, 63 of which fail on master (it was 90 and 52). The pairs are One correction to your pointer: the grid moved out of 4. — drafted by Claude Code on my behalf |
2 participants
Fixes #224.
fdistlogpdfwas evaluated aswhich breaks down at both ends: for a large
ν2the term(ν1 / ν2) * xis tiny and isdiscarded by the explicit
1 +, and for a largeν1the first and last terms both grow likeν1 * log(ν1)and cancel.It is now written in terms of the beta variate
u = ν1 * x / (ν1 * x + ν2), of which the densityis
u^a * (1 - u)^b / (x * beta(a, b))witha = ν1 / 2andb = ν2 / 2. Which expression isused depends on the ratio
r = u / (1 - u) = (ν1 / ν2) * x:1 / rfinite:log(u) = -log1p(1 / r)andlog(1 - u) = -log1p(r). This is symmetric in thetwo degrees of freedom, so neither failure mode can occur - a tiny
ν1 * x / ν2is not absorbedby an explicit
1 +, and no two terms growing likeν1 * log(ν1)have to cancel. The conditionis
1 / rrather than the more obviousr > 1because that is all the form actually requires,and the wider branch also covers the ridge
x ≈ ν2 / ν1, where the textbook form distributesaover two logarithms that cancel (1.3e-14 there, against 3.0e-16 here).ris zero or subnormal: the previous textbook expression, but withlog1p(r)inplace of
log(1 + r), which would drop a tinyrentirely. Splittinglog(u)intolog(ν1 / ν2) + log(x)rather than writinglog(r)keeps this branch accurate whenrunderflows, and it covers
x == 0as well asNaN, so no separate branch is needed for those.The ratio is formed as
(ν1 / ν2) * x, notν1 * x / ν2, so that an ordinaryxnearfloatmaxdoes not overflow it, and where
roverflows anywaylog1p(r)is taken aslog(ν1 / ν2) + log(x)- both summands are positive exactly when that happens.x < 0returns beforelog(x)is reached, which also removes themax(x, 0)that was onlythere to keep
logfrom erroring.Relative error of the density against a 512-bit reference:
The cdf family and the
ChainRulesCorerule are built on the same ratio. Both used to overflowat either end of the range, taking the result with them:
fdistlogcdf(3.0, 7.0, 1e-310)-Inf-1070.2806176752842fdistlogcdf(Float16(1), Float16(7), Float16(1e-4))-Inf-4.867fdistlogpdf(1.0, 0.5, 1e308)-Inf-888.325654114039R agrees to the last digit shown on both cdf values. The beta variate of the cdf is built as
r / (1 + r)rather thany / (y + ν2), whose denominator overflows forν2 > floatmax / 2.There the ratio is formed as
(ν1 * x) / ν2, the other way round fromfdistlogpdf: anoverflowing
ν1 * xonly saturatesuto 1, which is the right answer, and dividing afterwardskeeps
uexact whereν1 / ν2would go subnormal.The AD rule
∂/∂xwas((ν1 - 2) / x - ν1 * (ν1 + ν2) / (ν1 * x + ν2)) / 2, a difference of two terms bothgrowing like
ν1 / x, which cancels catastrophically. Combining them overν1 * x + ν2collapsesthe two into a single quotient, since
a * ν2 - b * ν1 * xis(ν1 * ν2 / 2) * (1 - x):The product is grouped as
ν1 * (ν2 * a)rather thanν1 * ν2 * a, becauseν1 * ν2overflows fordegrees of freedom of
floatmaxscale even where the derivative itself is ordinary, whereasν2 * ais a ratio of the two and stays in range. It is also the subexpression∂/∂ν1alreadycomputes.
0.0, true-1.1250.0, true-1.01.0e-308, true-5.1e-307NaN, true-5.0e99-Inf, true-1.25e199-Inf, true-8.3e198What kept this rewrite out until now is that it changes the value at
x = Inf. That turns out tobe a property of
a = (x - 1) / (ν1 * x + ν2), which isNaNthere and wrong once thedenominator overflows, rather than of the rewrite. With
aguarded, the gradient atx = Infgoes from
[NaN, NaN, 0.0]to[-0.4319, -Inf, -0.0], all three of which are the correct limits,and the other two partials become accurate for a large
xas well: at(3, 7, 1e308)thegradient goes from
[0.735, -Inf, 5.0e-309]to[-0.432, -353.5, -4.5e-308].Behaviour changes
fdistlogpdf(3, 7, Inf)returns-Infinstead ofNaN.x = Infis finite inν1and-Infinν2, where both wereNaN.NaNstill propagates, now throughxlogy.Tests
The
BigFloatreferences live intest/utils.jlnext toUtils; the assertions join the othernumerical issuetestsets intest/misc.jl, and the boundary values are a testset ofMiscalongside
support of binomial distribution. They check against aBigFloatevaluation of thetextbook formula rather than Rmath -
Rmath.dfis itself off by 2.5e-8 atν1 = 1e12, so itcannot pin that direction.
The references evaluate at 1024 bits rather than the default 256. That is not belt and braces: the
logbeta(a, b)of the textbook formula cancelslgamma(b) - lgamma(a + b)at the hundredth digitwhen the two degrees of freedom are 100 orders of magnitude apart, and
∂/∂xcancels to 300 digits.At 256 bits the reference for
logpdf(1e200, 1e300, 0.5)is off by a relative 3.6e3 and the one fordlogpdf_dx(1e300, 1e200, 2.0)by a relative 1.0 - either would have validated anything.The grid covers both large-degree-of-freedom directions over the whole range of
x, from wherethe ratio underflows to where it overflows, including the ridge
x ≈ ν2 / ν1, and both degrees offreedom at
floatmaxscale, where the productsν1 * xandν1 * ν2overflow. 63 of its 106assertions fail on master; the rest are the points master already gets right, kept so that the
grid is a grid. The ridge is not one of the 63 - master's error there is 1.3e-14, comfortably
inside the 1e-12 tolerance - so those points are coverage of the branch boundary rather than a
pin on its accuracy. Tightening the tolerance far enough to pin it would put the tolerance within
a factor of two of master's error, which is exactly the sort of margin that turns into a spurious
CI failure on another platform. For the same reason the
floatmax-scale pairs are1e200against1e300, accurate to 7.6e-14, and notfloatmaxitself against1e300, where thelogbetaof theimplementation loses enough to reach 4.9e-13. The cdf is pinned against the
u -> 0asymptoteu^a / (a * beta(a, b)), which needs no incomplete beta function and is exact to 2.1e-16 atx = 1e-310.The boundary testset covers the sign of
x,±Inf,NaN,x == 0and-0.0; a separatetestset checks that all three branches of
fdistlogpdfand the four cdf functions returnfloat(promote_type(...))for mixedInt,Rational,Float16,Float32andBigIntarguments. Moving the
x < 0result off the shared ternary and onto its own branch exposed thefact that it was never executed: there were no fdist tests below the support until #209 added
some, and the ternary counted as covered whenever
x > 0ran.For the AD rule,
test_fruleis used at two arguments where it used to overflow, with the tangentof
xscaled toxitself since the default tangent is of order one and its finite differenceswould leave the domain.
∂/∂xcannot be checked that way at all - a finite difference inν1islost in the rounding of
ν1itself - so the pullback is evaluated directly against aBigFloatreference.
#225 tracks the remaining case where both degrees of freedom are large, which this does not
address.
🤖 Generated with Claude Code