fix: stop destroying tail precision in evidence and truncation
`erfc` is sound — it holds ~1e-7 *relative* accuracy down to 1e-296 with no tail degradation. Three expressions built on it threw that away by subtracting quantities that both approach the same value. 1. `cavity_evidence` computed `1.0 - cdf(margin, ..)`, which is algebraically `sf(margin, ..)` and numerically a catastrophe: 7% error by eight sigma, and exactly zero past ~8.3, where the true probability is 1e-19 and perfectly representable. Clamped, that reached `log_evidence` as ln(f64::MIN_POSITIVE) = -708 whatever the truth was — off by 665 nats at nine sigma. `1 - cdf` is smallest precisely when the result contradicts the prior, so the model-comparison number was worst for upsets: the observation it exists to notice. Adds `sf`, the survival function, computed without the subtraction. The tie branch picks whichever tail keeps both of its terms small, for the same reason. 2. `v_w` computed the inverse Mills ratio as `pdf(-a) / cdf(-a)`. Both underflow together past about 39 sigma, giving `0 / 0` and putting NaN straight into the posterior. Adds `erfcx`, so the shared `exp(-alpha^2 / 2)` cancels analytically instead of being evaluated twice and divided. 3. With that fixed, `w = v * (v - alpha)` became the next casualty: `v` tends to `alpha`, so the gap lost every digit and drove `w` above 1, making `sqrt(1 - w)` NaN at alpha = 1e6. The gap now comes from its asymptotic series, which forms no difference at all. The tie branch had the same defect one expression over — `v * v - u` with both terms at 1e18 returned w = -128 — and a far-tail window is indistinguishable from a half-line, so it shares the asymptotic. No public signature changes, and no existing golden moved: every one of these only alters regions the old code got wrong. The two identity tests are asserted at 1e-6 rather than tighter because `erfc` is not exactly antisymmetric — `erfc(z) + erfc(-z)` differs from 2 by ~3e-8, and `erfc(0)` returns 1.00000003. That floor is tracked separately. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ
This commit is contained in:
Claude Opus 5
+67
-7
@@ -2,6 +2,7 @@ use crate::{
|
||||
N_INF, approx, cdf,
|
||||
factor::{Factor, VarId, VarStore},
|
||||
gaussian::Gaussian,
|
||||
sf,
|
||||
};
|
||||
|
||||
/// EP truncation factor on a diff variable.
|
||||
@@ -74,16 +75,29 @@ impl Factor for TruncFactor {
|
||||
|
||||
/// P(diff > margin) for non-tie, P(|diff| < margin) for tie.
|
||||
///
|
||||
/// Clamped to a positive floor: for a near-certain outcome the tail rounds to
|
||||
/// exactly 0.0, and the `erfc` approximation used by `cdf` carries ~1e-7 error
|
||||
/// so it can even return slightly more than 1.0, making the difference
|
||||
/// negative. Either would send `log_evidence` to `-inf` or NaN and poison the
|
||||
/// sum across the whole history.
|
||||
/// Both branches pick whichever tail keeps their terms *small*, because the
|
||||
/// alternative is subtracting two numbers that both approach 1. That
|
||||
/// subtraction is not a rounding detail: it loses every digit of an unlikely
|
||||
/// outcome's evidence, and an unlikely outcome is precisely the one worth
|
||||
/// scoring. `1 - cdf` returned exactly zero past ~8.3 sigma, where the true
|
||||
/// probability is 1e-19; clamped, that reached `log_evidence` as -708 instead
|
||||
/// of -43.
|
||||
///
|
||||
/// The clamp remains as a guard rather than a workaround: `erfc` carries ~1e-7
|
||||
/// relative error, so a probability of exactly 1 can still come back a hair
|
||||
/// above it, and `ln` of a negative would poison the sum for the whole history.
|
||||
fn cavity_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 {
|
||||
let (mu, sigma) = (diff.mu(), diff.sigma());
|
||||
|
||||
let raw = if tie {
|
||||
cdf(margin, diff.mu(), diff.sigma()) - cdf(-margin, diff.mu(), diff.sigma())
|
||||
if mu < -margin {
|
||||
// Both CDFs sit against 1 here; both survival terms are small.
|
||||
sf(-margin, mu, sigma) - sf(margin, mu, sigma)
|
||||
} else {
|
||||
1.0 - cdf(margin, diff.mu(), diff.sigma())
|
||||
cdf(margin, mu, sigma) - cdf(-margin, mu, sigma)
|
||||
}
|
||||
} else {
|
||||
sf(margin, mu, sigma)
|
||||
};
|
||||
|
||||
raw.clamp(f64::MIN_POSITIVE, 1.0)
|
||||
@@ -132,6 +146,52 @@ mod tests {
|
||||
assert_eq!(f.evidence_cached.unwrap(), first);
|
||||
}
|
||||
|
||||
/// The defect this guards: `1 - cdf` collapsed to zero for a surprising
|
||||
/// result, the clamp turned that into `f64::MIN_POSITIVE`, and
|
||||
/// `log_evidence` reported ln of *that* — about -708 whatever the truth
|
||||
/// was. An upset is the observation a model-comparison score exists to
|
||||
/// notice, so it was wrong exactly where it mattered.
|
||||
#[test]
|
||||
fn evidence_of_an_upset_is_not_flattened_to_the_clamp_floor() {
|
||||
// diff ~ N(-9, 1) with margin 0: the favoured side lost by nine sigma.
|
||||
let evidence = cavity_evidence(Gaussian::from_ms(-9.0, 1.0), 0.0, false);
|
||||
|
||||
assert!(
|
||||
evidence > f64::MIN_POSITIVE,
|
||||
"evidence collapsed onto the clamp floor: {evidence}"
|
||||
);
|
||||
// P(X > 0) for X ~ N(-9, 1) is the standard normal tail at 9 sigma.
|
||||
assert!(
|
||||
(evidence - 1.128_588e-19).abs() / 1.128_588e-19 < 1e-6,
|
||||
"expected ~1.13e-19, got {evidence}"
|
||||
);
|
||||
assert!(
|
||||
(evidence.ln() + 43.628).abs() < 1e-2,
|
||||
"log evidence {} should be about -43.6, not -708",
|
||||
evidence.ln()
|
||||
);
|
||||
}
|
||||
|
||||
/// Evidence must stay finite and positive however extreme the mismatch,
|
||||
/// since `log_evidence` sums across the whole history and one `-inf` or
|
||||
/// `NaN` poisons all of it.
|
||||
#[test]
|
||||
fn evidence_stays_positive_and_finite_at_any_separation() {
|
||||
for mu in [-300.0f64, -50.0, -9.0, 0.0, 9.0, 50.0, 300.0] {
|
||||
for tie in [false, true] {
|
||||
let e = cavity_evidence(Gaussian::from_ms(mu, 1.0), 1.0, tie);
|
||||
assert!(
|
||||
e.is_finite() && e > 0.0 && e <= 1.0,
|
||||
"mu={mu} tie={tie}: evidence {e} is not a probability"
|
||||
);
|
||||
assert!(
|
||||
e.ln().is_finite(),
|
||||
"mu={mu} tie={tie}: ln evidence is not finite"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tie_evidence_uses_two_sided() {
|
||||
let mut vars = VarStore::new();
|
||||
|
||||
+287
-8
@@ -168,6 +168,21 @@ pub const ITERATIONS: usize = 30;
|
||||
pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION;
|
||||
|
||||
const SQRT_TAU: f64 = 2.5066282746310002;
|
||||
/// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction.
|
||||
const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3;
|
||||
/// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form.
|
||||
const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4;
|
||||
/// How many window widths into the tail before a tie window is treated as a
|
||||
/// half-line. Beyond this the truncated mass is concentrated within `1/alpha`
|
||||
/// of the near edge, so the far edge contributes nothing measurable.
|
||||
const HALF_LINE_WINDOW: f64 = 10.0;
|
||||
/// Where `v - alpha` switches from subtraction to its asymptotic series.
|
||||
///
|
||||
/// The subtraction loses roughly `eps * alpha^2` of relative precision, and the
|
||||
/// four-term series is good to ~1e-10 by here, so the two are at their closest
|
||||
/// agreement around this point. Below it the subtraction is exact; above it the
|
||||
/// series is.
|
||||
const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;
|
||||
|
||||
pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);
|
||||
pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0);
|
||||
@@ -259,6 +274,50 @@ pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
|
||||
0.5 * erfc(z)
|
||||
}
|
||||
|
||||
/// `P(X > x)` for `X ~ N(mu, sigma^2)`.
|
||||
///
|
||||
/// The survival function, computed directly rather than as `1 - cdf(..)`.
|
||||
///
|
||||
/// The two are algebraically identical and numerically are not. `cdf` returns
|
||||
/// a value approaching 1 for an upper tail, so subtracting it from 1 cancels
|
||||
/// away every significant digit the tail had: measured against this function,
|
||||
/// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly
|
||||
/// zero beyond about 8.3 sigma — where the true value is still 1e-19 and
|
||||
/// perfectly representable. `erfc` itself holds ~1e-7 *relative* accuracy down
|
||||
/// to 1e-296, so the precision is there to keep; only the subtraction threw it
|
||||
/// away.
|
||||
///
|
||||
/// This matters most where evidence is smallest, which is exactly where an
|
||||
/// upset makes it interesting: `ln` of a clamped zero is -708 regardless of
|
||||
/// whether the truth was -43 or -600.
|
||||
pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 {
|
||||
0.5 * erfc((x - mu) / (sigma * SQRT_2))
|
||||
}
|
||||
|
||||
/// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`.
|
||||
///
|
||||
/// Exists so the exponential factor common to a Gaussian density and its tail
|
||||
/// integral can be cancelled *analytically* instead of being computed twice
|
||||
/// and divided. Both underflow to zero past about 26 sigma, and their ratio is
|
||||
/// then `0/0` — finite in the limit, `NaN` in floating point.
|
||||
fn erfcx(x: f64) -> f64 {
|
||||
if x < 2.0 {
|
||||
// Below the crossover neither factor is extreme: erfc is O(1) and
|
||||
// exp(x^2) is at most e^4, so the direct product is exact enough and
|
||||
// cheaper than the continued fraction.
|
||||
(x * x).exp() * erfc(x)
|
||||
} else {
|
||||
// erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))),
|
||||
// evaluated by backward recurrence. Converges quickly for x >= 2 and,
|
||||
// unlike the product form, never touches an exponential.
|
||||
let mut f = 0.0;
|
||||
for n in (1..=60u32).rev() {
|
||||
f = (f64::from(n) * 0.5) / (x + f);
|
||||
}
|
||||
FRAC_1_SQRT_PI / (x + f)
|
||||
}
|
||||
}
|
||||
|
||||
fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
|
||||
let normalizer = (SQRT_TAU * sigma).powi(-1);
|
||||
let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp();
|
||||
@@ -266,25 +325,100 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
|
||||
normalizer * functional
|
||||
}
|
||||
|
||||
/// Truncated-Gaussian correction terms `(v, w)`.
|
||||
///
|
||||
/// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose
|
||||
/// numerator and denominator underflow together in the tails, so both are
|
||||
/// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled
|
||||
/// analytically rather than evaluated and divided out. Without that, a
|
||||
/// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN`
|
||||
/// straight into the posterior.
|
||||
/// Truncation terms for a boundary `alpha` standard deviations into the upper
|
||||
/// tail, from the asymptotic expansion of the inverse Mills ratio.
|
||||
///
|
||||
/// `v` tends to `alpha` out here, so the gap between them cannot be obtained by
|
||||
/// subtracting one from the other — the series computes the gap directly, and
|
||||
/// `w = v * gap` then never forms the difference of two large near-equal
|
||||
/// numbers. A far-tail *window* behaves like a half-line once it is more than a
|
||||
/// few multiples of its own width from the mean, so the tie branch shares this.
|
||||
fn half_line_truncation(alpha: f64) -> (f64, f64) {
|
||||
let inv = alpha.recip();
|
||||
let inv_sq = inv * inv;
|
||||
let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
|
||||
let v = alpha + gap;
|
||||
|
||||
(v, v * gap)
|
||||
}
|
||||
|
||||
fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
|
||||
if !tie {
|
||||
let alpha = (margin - mu) / sigma;
|
||||
|
||||
let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
|
||||
let w = v * (v + (-alpha));
|
||||
// v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs
|
||||
// the gap `v - alpha` as well as v itself. Far into the tail v tends to
|
||||
// alpha, so that gap is a subtraction of two nearly equal numbers and
|
||||
// loses every digit it has: at alpha = 1e6 it drove w above 1 and made
|
||||
// `sqrt(1 - w)` NaN. Past the crossover the gap comes from its
|
||||
// asymptotic series instead, which has no subtraction in it.
|
||||
if alpha >= ASYMPTOTIC_MILLS_ALPHA {
|
||||
return half_line_truncation(alpha);
|
||||
}
|
||||
|
||||
(v, w)
|
||||
let (v, gap) = if alpha > 0.0 {
|
||||
// Both terms carry exp(-alpha^2 / 2); in scaled form it cancels
|
||||
// and the result stays exact however far into the tail alpha sits.
|
||||
let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2);
|
||||
(v, v - alpha)
|
||||
} else {
|
||||
// Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing.
|
||||
let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
|
||||
(v, v - alpha)
|
||||
};
|
||||
|
||||
(v, v * gap)
|
||||
} else {
|
||||
// v is odd in mu and w is even, so fold to mu <= 0. Both truncation
|
||||
// points then sit in the upper tail, where the scaled form applies.
|
||||
let flipped = mu > 0.0;
|
||||
let mu = if flipped { -mu } else { mu };
|
||||
|
||||
let alpha = (-margin - mu) / sigma;
|
||||
let beta = (margin - mu) / sigma;
|
||||
|
||||
let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0))
|
||||
/ (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
|
||||
let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0))
|
||||
/ (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
|
||||
// `w` comes out of `v * v - u`, and both terms grow as alpha^2 while
|
||||
// their difference stays O(1) — at alpha = 1e9 that subtraction had no
|
||||
// digits left and returned w = -128, making `sqrt(1 - w)` nonsense.
|
||||
// Once the window sits many of its own widths into the tail it is
|
||||
// indistinguishable from a half-line, so the asymptotic covers it with
|
||||
// no subtraction at all.
|
||||
if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW {
|
||||
let (v, w) = half_line_truncation(alpha);
|
||||
return (if flipped { -v } else { v }, w);
|
||||
}
|
||||
|
||||
let (v, u) = if alpha > 0.0 {
|
||||
// beta > alpha > 0, so this ratio of exponentials is at most 1 and
|
||||
// cannot overflow.
|
||||
let scale = (0.5 * (alpha * alpha - beta * beta)).exp();
|
||||
let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));
|
||||
|
||||
(
|
||||
(1.0 - scale) / SQRT_TAU / denominator,
|
||||
(alpha - beta * scale) / SQRT_TAU / denominator,
|
||||
)
|
||||
} else {
|
||||
// The interval straddles the mean, so nothing here is small.
|
||||
let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0);
|
||||
|
||||
(
|
||||
(pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator,
|
||||
(alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator,
|
||||
)
|
||||
};
|
||||
|
||||
let w = -(u - v.powi(2));
|
||||
|
||||
(v, w)
|
||||
(if flipped { -v } else { v }, w)
|
||||
}
|
||||
}
|
||||
|
||||
@@ -467,6 +601,151 @@ mod tests {
|
||||
assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]);
|
||||
}
|
||||
|
||||
/// Upper-tail values of the standard normal, from published tables. The
|
||||
/// point is not the digits — `erfc` only carries ~1e-7 relative — but that
|
||||
/// a number comes back at all: `1 - cdf` returned exactly zero for every
|
||||
/// one of these.
|
||||
#[test]
|
||||
fn survival_function_survives_the_far_tail() {
|
||||
for (z, expected) in [
|
||||
(9.0f64, 1.128_588e-19),
|
||||
(12.0, 1.776_482e-33),
|
||||
(20.0, 2.753_624e-89),
|
||||
(37.0, 5.725_571e-300),
|
||||
] {
|
||||
let got = sf(z, 0.0, 1.0);
|
||||
assert!(got > 0.0, "sf({z}) collapsed to zero");
|
||||
assert!(
|
||||
(got - expected).abs() / expected < 1e-6,
|
||||
"sf({z}) = {got}, expected ~{expected}"
|
||||
);
|
||||
assert_eq!(
|
||||
1.0 - cdf(z, 0.0, 1.0),
|
||||
0.0,
|
||||
"the naive form should still be zero here"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Where no cancellation happens the two forms must agree exactly enough
|
||||
/// that nothing else in the crate shifts.
|
||||
#[test]
|
||||
fn survival_function_matches_the_naive_form_where_that_form_works() {
|
||||
for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] {
|
||||
let naive = 1.0 - cdf(z, 0.0, 1.0);
|
||||
let direct = sf(z, 0.0, 1.0);
|
||||
// Bounded by `erfc`'s own ~1e-7 relative error, not by the
|
||||
// subtraction: the two forms evaluate `erfc` at different points
|
||||
// and the approximation is not exactly antisymmetric.
|
||||
assert!(
|
||||
(naive - direct).abs() < 1e-6,
|
||||
"z={z}: naive {naive} vs direct {direct}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn survival_and_cdf_partition_the_mass() {
|
||||
for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] {
|
||||
let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0);
|
||||
// `erfc(z) + erfc(-z) == 2` only to the accuracy of the
|
||||
// approximation, which is ~1e-7 relative.
|
||||
assert!((total - 1.0).abs() < 1e-6, "z={z}: {total}");
|
||||
}
|
||||
}
|
||||
|
||||
/// `erfcx` switches formulation at x = 2; the two sides must meet.
|
||||
#[test]
|
||||
fn erfcx_is_continuous_across_its_crossover() {
|
||||
for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] {
|
||||
let direct = (x * x).exp() * erfc(x);
|
||||
let scaled = erfcx(x);
|
||||
assert!(
|
||||
(direct - scaled).abs() / scaled < 1e-6,
|
||||
"x={x}: direct {direct} vs erfcx {scaled}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly
|
||||
/// where `exp(x^2)` overflows and `erfc(x)` underflows.
|
||||
#[test]
|
||||
fn erfcx_stays_finite_where_its_factors_do_not() {
|
||||
for x in [27.0f64, 50.0, 1.0e3, 1.0e8] {
|
||||
let scaled = erfcx(x);
|
||||
assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}");
|
||||
// Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)).
|
||||
let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt());
|
||||
assert!(
|
||||
(scaled - asymptote).abs() / asymptote < 1e-2,
|
||||
"erfcx({x}) = {scaled} strays from its asymptote {asymptote}"
|
||||
);
|
||||
assert!(
|
||||
(x * x).exp().is_infinite(),
|
||||
"x={x} should overflow the direct form"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Truncation must never produce a non-finite posterior. Before the scaled
|
||||
/// formulation these returned NaN from `0 / 0` past about 39 sigma.
|
||||
#[test]
|
||||
fn truncation_stays_finite_arbitrarily_far_into_the_tail() {
|
||||
for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] {
|
||||
for tie in [false, true] {
|
||||
let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
|
||||
assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}");
|
||||
assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}");
|
||||
// sigma_trunc = sigma * sqrt(1 - w) must stay real.
|
||||
assert!(
|
||||
(0.0..=1.0).contains(&w),
|
||||
"alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary"
|
||||
);
|
||||
|
||||
let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
|
||||
assert!(
|
||||
mu_t.is_finite() && sigma_t.is_finite(),
|
||||
"alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The Mills gap switches from subtraction to series at alpha = 100. Both
|
||||
/// are supposed to be right there; if they disagree, the crossover is in
|
||||
/// the wrong place.
|
||||
#[test]
|
||||
fn the_mills_gap_series_meets_the_scaled_form() {
|
||||
for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] {
|
||||
let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha;
|
||||
let inv = alpha.recip();
|
||||
let inv_sq = inv * inv;
|
||||
let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
|
||||
assert!(
|
||||
(scaled - series).abs() / series < 1e-9,
|
||||
"alpha={alpha}: scaled {scaled} vs series {series}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and
|
||||
/// w is even. Assert the symmetry the implementation relies on.
|
||||
#[test]
|
||||
fn tie_truncation_is_odd_in_v_and_even_in_w() {
|
||||
for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] {
|
||||
let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true);
|
||||
let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true);
|
||||
assert!(
|
||||
(v_pos + v_neg).abs() < 1e-9,
|
||||
"mu={mu}: v should be odd, got {v_pos} and {v_neg}"
|
||||
);
|
||||
assert!(
|
||||
(w_pos - w_neg).abs() < 1e-9,
|
||||
"mu={mu}: w should be even, got {w_pos} and {w_neg}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_quality() {
|
||||
let a = Gaussian::from_ms(25.0, 3.0);
|
||||
|
||||
Reference in New Issue
Block a user