back to table · edit · history · where entries came from · files · download
5674 bytes, as of the version from 2026-10-16 06:19 (current). Recorded here, not run.
"""Quantiles of Student's t-distribution -- numberdb.org/T454.
Run with SageMath (or compatible passagemath):
sage -pip install numberdb
sage -python generate.py # verify all stored entries; key for draft
sage -python generate.py --fill # write entries to the existing draft
NUMBERDB_API_KEY authorizes access to the draft. --fill writes numbers and
attachments; it does not publish the table or offer it for review.
The central, location-zero, scale-one distribution is used. Parameters are
(nu, p), with exact rational lower-tail probability p. The 14 probabilities
are the family's conventional significance levels and quartiles. Degrees of
freedom 1..30 cover small samples; 40, 50, 60, 100 are larger reference rows,
the same selection as the chi-squared companion T453.
Each positive quantile is enclosed by rational bisection. The CDF is evaluated
in Arb via I_z(a,b) = z^a 2F1(a,1-b;a+1;z)/(a B(a,b)), DLMF 8.17.7.
The beta complement keeps the hypergeometric argument at most 1/2. Strict
ball comparisons certify every bracket update. The returned ball is the union
of both endpoints, never a rounded midpoint. Symmetry supplies negative rows.
All medians and the Cauchy quartiles +/-1 are returned exactly.
WORKING_GUARD=160 leaves ample space for subtraction and beta evaluation;
checks.py measures the minimum relative accuracy of every returned enclosure.
The stopping bracket has at least 12 guard decimal places beyond the 100
requested digits. Independent checks integrate cos(theta)^(nu-1) by a finite
recurrence, sharing neither beta nor hypergeometric code with this generator.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
WORKING_GUARD = 160
DEGREES = tuple(range(1, 31)) + (40, 50, 60, 100)
PROBABILITIES = tuple(QQ(x) for x in (
'1/200', '1/100', '1/40', '1/20', '1/10', '1/4', '1/2',
'3/4', '9/10', '19/20', '39/40', '99/100', '199/200', '999/1000'))
@lru_cache(maxsize=None)
def beta_normalizer(a, b, prec):
C = ComplexBallField(prec)
return C(a).gamma() * C(b).gamma() / C(a + b).gamma()
def beta_regularized(z, a, b, prec):
C = ComplexBallField(prec)
z, aa, bb = C(z), C(a), C(b)
ans = z**aa * z.hypergeometric([aa, 1-bb], [aa+1])
ans /= aa * beta_normalizer(a, b, prec)
assert abs(ans.real()) < 2 and ans.imag().contains_zero()
return ans.real()
def cdf_positive(nu, x, prec):
"""Certified CDF at a nonnegative exact rational x."""
R = RealBallField(prec)
x = QQ(x)
assert x >= 0
if x == 0:
return R(QQ(1)/2)
z = QQ(nu) / (nu + x*x)
if z <= QQ(1)/2:
return 1 - beta_regularized(z, QQ(nu)/2, QQ(1)/2, prec)/2
return (1 + beta_regularized(1-z, QQ(1)/2, QQ(nu)/2, prec))/2
@lru_cache(maxsize=None)
def positive_quantile(nu, p, digits):
prec = numberdb.bits(digits, losing=WORKING_GUARD)
R = RealBallField(prec)
lo, hi = QQ(0), QQ(1)
while not cdf_positive(nu, hi, prec) > R(p):
hi *= 2
if hi > 2**30:
raise ArithmeticError('Could not certify an upper bracket')
tolerance = QQ(1) / (10**(digits + 12))
for _ in range(4*prec):
if hi-lo < tolerance * max(QQ(1), lo):
value = R(lo).union(R(hi))
assert 0 < value < 2**30
assert R(value.rad()) < R(abs(value.center())) * R(QQ(1)/10**(digits+10))
return value
mid = (lo+hi)/2
fm = cdf_positive(nu, mid, prec)
if fm < R(p):
lo = mid
elif fm > R(p):
hi = mid
else:
# Never select a side from an overlapping CDF ball.
# Two quarter points preserve a bracket even if the midpoint
# itself is the exact root.
left, right = (lo+mid)/2, (mid+hi)/2
fl, fr = cdf_positive(nu, left, prec), cdf_positive(nu, right, prec)
if fl < R(p) and fr > R(p):
lo, hi = left, right
else:
raise ArithmeticError('CDF uncertainty exceeds bracket resolution')
raise ArithmeticError('Bisection failed to reach the requested precision')
class StudentQuantiles(numberdb.Generator):
table = 'T454'
parameters = ('nu', 'p')
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'nist-critical-values.json', 'check.py',
'extract_nist.py', 'nist-source.html')
def enumerate(self):
for nu in DEGREES:
for p in PROBABILITIES:
yield {'nu': str(nu), 'p': str(p)}
def value(self, params, digits):
nu, p = int(params['nu']), QQ(params['p'])
assert nu > 0 and 0 < p < 1
if p == QQ(1)/2:
return ZZ(0)
if nu == 1 and p in (QQ(1)/4, QQ(3)/4):
return ZZ(1 if p > QQ(1)/2 else -1)
if p < QQ(1)/2:
return -positive_quantile(nu, 1-p, digits)
return positive_quantile(nu, p, digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator = StudentQuantiles()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill certified central Student t lower-tail quantiles',
assisted_by='codex-cli'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)