back to table · edit · history · where entries came from · files · download
7184 bytes, as of the version from 2026-10-16 06:19 (current). Recorded here, not run.
"""Independent checks for Student t quantiles -- numberdb.org/T454.
sage -python checks.py # check locally generated values
sage -python checks.py --stored # read the draft; needs NUMBERDB_API_KEY
The independent CDF uses x=sqrt(nu)*tan(theta), reducing the density integral
to J_m(theta)=integral_0^theta cos(u)^m du, m=nu-1. Integration by parts gives
J_m=sin(theta)cos(theta)^(m-1)/m+(m-1)J_(m-2)/m. Base cases are J_0=theta
and J_1=sin(theta). Its normalization obeys c_1=1/pi, c_2=1/2,
c_(nu+2)=(nu+1)c_nu/nu. No hypergeometric, beta or gamma evaluation is used.
Each stored decimal denotes plus or minus one last-place unit. At both
endpoints we prove the independent CDF is on the correct side of p. Finite,
narrow ball checks precede every numerical comparison. Exact rows are
validated separately from mathematical symmetry and Cauchy quartiles.
The NIST comparison tests all six shared probabilities in all 34 degree rows.
"""
import json
import os
import sys
from pathlib import Path
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from numberdb._write import to_text
from check import exactness, measure, names_its_rings
from generate import StudentQuantiles, DEGREES, PROBABILITIES
R = RealBallField(600)
def rational_decimal(s):
f = Fraction(str(s))
return QQ(f.numerator)/QQ(f.denominator)
def endpoints(s):
s = str(s)
q = rational_decimal(s)
if '.' not in s and 'e' not in s.lower():
return q, q
mantissa, _, exponent = s.lower().partition('e')
places = len(mantissa.split('.')[1]) if '.' in mantissa else 0
ulp = QQ(10)**(int(exponent or '0')-places)
return q-ulp, q+ulp
def finite_narrow(v, tolerance=QQ(1)/10**130):
assert v.is_finite() and v.rad() < tolerance, str(v)
return v
def independent_cdf(nu, x):
x = R(x)
den = (R(nu)+x*x).sqrt()
sin_theta, cos_theta = x/den, R(nu).sqrt()/den
if nu % 2:
j = (x/R(nu).sqrt()).arctan()
c = 1/R.pi()
start = 2
else:
j = sin_theta
c = R(QQ(1)/2)
start = 3
for m in range(start, nu, 2):
j = (sin_theta*cos_theta**(m-1)+(m-1)*j)/m
for k in range(1 if nu % 2 else 2, nu, 2):
c *= R(QQ(k+1)/k)
return R(QQ(1)/2)+c*j
def controls():
for nu in DEGREES:
f = finite_narrow(independent_cdf(nu, 0))
assert (f-R(QQ(1)/2)).contains_zero()
f = finite_narrow(independent_cdf(1, 1))
assert (f-R(QQ(3)/4)).contains_zero()
f = finite_narrow(independent_cdf(2, R(2).sqrt()))
expected = (1+1/R(2).sqrt())/2
assert finite_narrow(f-expected).contains_zero()
# A nontrivial positive control rejects a routine returning zero always.
assert f > R(QQ(4)/5)
def check_rows(rows):
controls()
wanted = {(str(n),str(p)) for n in DEGREES for p in PROBABILITIES}
assert set(rows) == wanted and len(rows) == 476
balls = {}
counts = {'entries':len(rows),'certified_decimal_brackets':0,'exact_rows':0,
'symmetry_pairs':0,'cauchy_closed_forms':0,'nu2_closed_forms':0,
'nist_values':0,'probability_order_pairs':0,'degree_order_pairs':0,
'two_sided_tail_checks':0}
max_nist_error = QQ(0)
nist = json.loads(Path(__file__).with_name('nist-critical-values.json').read_text())['values']
for (ns,ps),text in rows.items():
nu, p = int(ns), QQ(ps)
lo, hi = endpoints(text)
b = R(lo).union(R(hi))
assert b.is_finite() and b.rad() < QQ(1)/10**95
balls[ns,ps] = b
if lo == hi:
assert (p == QQ(1)/2 and lo == 0) or (nu == 1 and (p,lo) in ((QQ(1)/4,-1),(QQ(3)/4,1)))
assert finite_narrow(independent_cdf(nu,lo)-R(p)).contains_zero()
counts['exact_rows'] += 1
else:
left = finite_narrow(independent_cdf(nu,lo))
right = finite_narrow(independent_cdf(nu,hi))
assert left < R(p) < right, (ns,ps,text,str(left-R(p)),str(right-R(p)))
counts['certified_decimal_brackets'] += 1
if p > QQ(1)/2:
# P(-t <= T <= t)=1-alpha, alpha=2(1-p).
for x in (lo,hi):
fc = finite_narrow(independent_cdf(nu,x))
fm = finite_narrow(independent_cdf(nu,-x))
assert finite_narrow(fc+fm-1).contains_zero()
counts['two_sided_tail_checks'] += 1
if nu == 1:
y = -(R.pi()*R(p)).cos()/(R.pi()*R(p)).sin()
finite_narrow(y)
assert b.overlaps(y)
counts['cauchy_closed_forms'] += 1
if nu == 2:
y = R(2*p-1)/(2*R(p)*R(1-p)).sqrt()
finite_narrow(y)
assert b.overlaps(y)
counts['nu2_closed_forms'] += 1
if ps in nist[ns]:
error=abs(rational_decimal(text)-rational_decimal(nist[ns][ps]))
assert error <= QQ(1)/2000, (ns,ps,str(error))
max_nist_error=max(max_nist_error,error)
counts['nist_values'] += 1
for nu in DEGREES:
for p,q in zip(PROBABILITIES,PROBABILITIES[1:]):
assert balls[str(nu),str(p)] < balls[str(nu),str(q)]
counts['probability_order_pairs'] += 1
for p in PROBABILITIES:
if p < QQ(1)/2 and 1-p in PROBABILITIES:
assert balls[str(nu),str(p)].overlaps(-balls[str(nu),str(1-p)])
counts['symmetry_pairs'] += 1
for p in PROBABILITIES:
if p == QQ(1)/2: continue
for n,m in zip(DEGREES,DEGREES[1:]):
a,b=balls[str(n),str(p)],balls[str(m),str(p)]
assert a < b if p < QQ(1)/2 else a > b
counts['degree_order_pairs'] += 1
counts['largest_nist_absolute_error'] = str(float(max_nist_error))
counts['written_size'] = measure({','.join(k):v for k,v in rows.items()})
print('CHECKS_JSON '+json.dumps(counts,sort_keys=True))
return counts
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
if '--stored' in sys.argv or os.environ.get('CHECK_STORED') == '1':
doc=numberdb.table('T454')
numbers=doc['Numbers']
rows={(str(n),str(p)):str(v['number'] if isinstance(v,dict) else v)
for n,group in numbers.items() for p,v in group.items()}
else:
gen=StudentQuantiles(); values={}; rows={}; accuracies=[]
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
for params in gen.enumerate():
value=gen.value(params,100)
key=(params['nu'],params['p']); values[key]=value
rows[key]=to_text(value,100)
if hasattr(value,'rad'):
accuracy=-(R(value.rad())/abs(R(value.center()))).log()/R(2).log()
accuracies.append(float(accuracy.center()))
assert not exactness(values)
print('MINIMUM_RELATIVE_ACCURACY_BITS',min(accuracies))
check_rows(rows)
print('ENTRIES_JSON '+json.dumps([{'params':{'nu':n,'p':p},'number':v} for (n,p),v in rows.items()]))
if __name__ == '__main__':
main()