back to table · edit · history · where entries came from · files · download
5386 bytes, as of the version from 2026-10-16 06:58 (current). Recorded here, not run.
"""Quantiles of the F-distribution -- numberdb.org/T455.
Run with SageMath (or passagemath with the required ball-arithmetic modules):
sage -pip install numberdb
sage -python generate.py # verify every stored entry
sage -python generate.py --fill # fill the existing draft
For an unpublished draft, set NUMBERDB_API_KEY or pipe the key on stdin with
NUMBERDB_KEY_FROM_STDIN=1. --fill writes entries; it never publishes a table
or offers it for review. NUMBERDB_PUBLISH=1 is the batch-run equivalent.
Convention: the central F law is (U/d1)/(V/d2) for independent chi-squared
U,V. p is a lower-tail probability, written as an exact rational. The grid
contains d1,d2=1,...,10, the median, and the 1%, 2.5%, and 5% tails. These
are small-sample variance-ratio and ANOVA critical values, not a decimal grid.
Each nonrational value is enclosed by bisection of the beta CDF on [0,1],
using the hypergeometric formula DLMF 8.17.7 and complement symmetry. Only
certified signs change a bracket. Exact rational endpoints are transformed
by x=d2*z/(d1*(1-z)); the returned ball encloses the WHOLE bracket. Working
precision is bits(digits)+128; termination requires a relative bracket width
at most 2**(-bits(digits)-24). Checks measure the retained precision over the
whole grid. Exact roots in the d1=2 or d2=2 formulas, and equal-degree medians,
are returned as rationals, never rounded decimals.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfi import RealIntervalField
WORKING_GUARD = 128
BRACKET_GUARD = 24
PROBABILITIES = ('1/100', '1/40', '1/20', '1/2', '19/20', '39/40', '99/100')
DEGREES = tuple(range(1, 11))
def rational_power(base, exponent):
"""Return a rational power only if an exact integer-root test proves it."""
n = int(exponent.denominator())
try:
root = QQ(base.numerator().nth_root(n)) / base.denominator().nth_root(n)
except ValueError:
return None
return root ** exponent.numerator()
def exact_quantile(d1, d2, p):
if d1 == d2 and p == QQ(1)/2:
return QQ(1)
if d2 == 2:
z = rational_power(p, QQ(2)/d1)
if z is not None:
return QQ(2)/d1*z/(1-z)
if d1 == 2:
w = rational_power(1-p, -QQ(2)/d2)
if w is not None:
return QQ(d2)/2*(w-1)
return None
@lru_cache(None)
def beta_setup(d1, d2, prec):
C = ComplexBallField(prec)
a, b = C(QQ(d1)/2), C(QQ(d2)/2)
return C, a, b, (a+b).gamma()/(a.gamma()*b.gamma())
def beta_cdf(z, d1, d2, prec):
R = RealBallField(prec)
if z == 0 or z == 1:
return R(z)
if z > QQ(1)/2:
return 1-beta_cdf(1-z, d2, d1, prec)
C, a, b, norm = beta_setup(d1, d2, prec)
t = C(z)
value = t**a*t.hypergeometric([a, 1-b], [a+1])*norm/a
if not value.real().is_finite() or not value.imag().is_finite() or not value.imag().contains_zero():
raise ArithmeticError('Invalid incomplete-beta enclosure')
return value.real()
def quantile(d1, d2, p, digits=100):
d1, d2, p = int(d1), int(d2), QQ(p)
if min(d1, d2) < 1 or not 0 < p < 1:
raise ValueError('Positive integer degrees and 0<p<1 required')
exact = exact_quantile(d1, d2, p)
if exact is not None:
return exact
target = numberdb.bits(digits, losing=BRACKET_GUARD)
prec = numberdb.bits(digits, losing=WORKING_GUARD)
lo, hi = QQ(0), QQ(1)
for _ in range(4*target):
if lo > 0 and hi < 1:
xlo = QQ(d2)*lo/(d1*(1-lo))
xhi = QQ(d2)*hi/(d1*(1-hi))
if xhi-xlo <= xlo/2**target:
result = RealBallField(prec)(RealIntervalField(prec)(xlo, xhi))
assert result.is_finite() and result > 0
return result
mid = (lo+hi)/2
residual = beta_cdf(mid, d1, d2, prec)-RealBallField(prec)(p)
if residual < 0:
lo = mid
elif residual > 0:
hi = mid
else:
left, right = (lo+mid)/2, (mid+hi)/2
if beta_cdf(left, d1, d2, prec) < RealBallField(prec)(p) \
and beta_cdf(right, d1, d2, prec) > RealBallField(prec)(p):
lo, hi = left, right
else:
prec *= 2
raise ArithmeticError('Certified bisection did not converge')
class FQuantiles(numberdb.Generator):
table = 'T455'
parameters = ('d1', 'd2', 'p')
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
for d1 in DEGREES:
for d2 in DEGREES:
for p in PROBABILITIES:
yield {'d1': d1, 'd2': d2, 'p': p}
def value(self, params, digits):
return quantile(params['d1'], params['d2'], params['p'], digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator = FQuantiles()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill 700 certified lower-tail F quantiles'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)