back to table · edit · history · where entries came from · files · download
9149 bytes, as of the version from 2026-10-16 06:37. Recorded here, not run.
"""Independent checks for Quantiles of the F-distribution -- numberdb.org/T455.
Run: sage -python checks.py [--stored]. Install numberdb and scipy first.
The default checks freshly generated entries; --stored reads the table back.
A private draft needs NUMBERDB_API_KEY, or NUMBERDB_KEY_FROM_STDIN=1 and stdin.
Keep generate.py beside this file. The repository helper check.py is attached
so its exactness, measurement and named-import checks remain reproducible.
The independent CDF uses elementary half-integer beta bases and finite
integration-by-parts recurrences. It calls no hypergeometric or gamma function.
Every written decimal is interpreted as its claimed last-place interval;
strict opposite CDF signs at its endpoints independently certify the root.
Controls precede measurements, including known uniform and arcsine beta laws.
"""
import json
import os
import sys
import time
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from generate import FQuantiles, PROBABILITIES
from check import exactness, names_its_rings, measure, agrees_with
R = RealBallField(700)
def elementary_cdf(x, d1, d2):
"""The F CDF, independently evaluated without beta/gamma/hypergeometric."""
z = R(d1)*R(x)/(R(d1)*R(x)+d2)
a = QQ(1)/2 if d1 % 2 else QQ(1)
b = QQ(1)/2 if d2 % 2 else QQ(1)
if a == b == 1:
val, beta = z, R(1)
elif a == b == QQ(1)/2:
val, beta = 2*z.sqrt().arcsin()/R.pi(), R.pi()
elif a == QQ(1)/2:
val, beta = z.sqrt(), R(2)
else:
val, beta = 1-(1-z).sqrt(), R(2)
while a < QQ(d1)/2:
val -= z**a*(1-z)**b/(a*beta)
beta *= R(a)/(a+b)
a += 1
while b < QQ(d2)/2:
val += z**a*(1-z)**b/(b*beta)
beta *= R(b)/(a+b)
b += 1
assert val.is_finite()
return val
def endpoints(text):
if '.' not in text and 'e' not in text.lower():
q = QQ(text)
return q, q
d = Decimal(text)
q = QQ(Fraction(d))
ulp = QQ(10)**d.as_tuple().exponent
return q-ulp, q+ulp
def narrow_zero(x, tolerance=180):
assert x.is_finite() and x.contains_zero()
assert x.abs() < R(10)**(-tolerance)
def controls():
for p in (QQ(1)/100, QQ(1)/3, QQ(1)/2, QQ(19)/20):
narrow_zero(elementary_cdf(p/(1-p),2,2)-R(p))
narrow_zero(elementary_cdf(QQ(1)/3,1,1)-R(QQ(1)/3))
narrow_zero(elementary_cdf(3,1,1)-R(QQ(2)/3))
narrow_zero(elementary_cdf(199,2,2)-R(QQ(199)/200))
assert elementary_cdf(18,2,2) < R(QQ(19)/20)
assert elementary_cdf(20,2,2) > R(QQ(19)/20)
print('controls: uniform beta, arcsine beta, exact F(2,2), and nonzero signs pass',flush=True)
# NIST/SEMATECH e-Handbook, 1.3.6.7.3, first ten rows and columns,
# upper-tail alpha=0.05. Rows index denominator degrees; columns numerator.
NIST_95 = '''
161.448 199.500 215.707 224.583 230.162 233.986 236.768 238.882 240.543 241.882
18.513 19.000 19.164 19.247 19.296 19.330 19.353 19.371 19.385 19.396
10.128 9.552 9.277 9.117 9.013 8.941 8.887 8.845 8.812 8.786
7.709 6.944 6.591 6.388 6.256 6.163 6.094 6.041 5.999 5.964
6.608 5.786 5.409 5.192 5.050 4.950 4.876 4.818 4.772 4.735
5.987 5.143 4.757 4.534 4.387 4.284 4.207 4.147 4.099 4.060
5.591 4.737 4.347 4.120 3.972 3.866 3.787 3.726 3.677 3.637
5.318 4.459 4.066 3.838 3.687 3.581 3.500 3.438 3.388 3.347
5.117 4.256 3.863 3.633 3.482 3.374 3.293 3.230 3.179 3.137
4.965 4.103 3.708 3.478 3.326 3.217 3.135 3.072 3.020 2.978
'''
def inspect_entries(records):
from scipy.stats import f, t
controls()
values, bounds = {}, {}
exact_count = 0
max_scipy = 0.0
for entry in records:
params = entry['params']
d1, d2, p = int(params['d1']), int(params['d2']), QQ(params['p'])
key = (d1, d2, str(p))
lo, hi = endpoints(entry['number'])
assert lo > 0 and hi >= lo
bounds[key] = (lo, hi)
values[key] = R(RealIntervalField(700)(lo,hi))
left, right = elementary_cdf(lo,d1,d2)-R(p), elementary_cdf(hi,d1,d2)-R(p)
if lo == hi:
exact_count += 1
narrow_zero(left)
else:
assert left < 0 < right, (key,'CDF does not bracket p')
assert (hi-lo)/lo < QQ(1)/10**98
estimate = f.ppf(float(p),d1,d2)
relerr = abs(estimate-float((lo+hi)/2))/float((lo+hi)/2)
max_scipy = max(max_scipy,relerr)
assert relerr < 3e-12, (key,relerr)
assert len(bounds)==700
assert len(bounds)==len(records)
print('independent endpoint checks:',len(records)-exact_count,'decimal enclosures;',exact_count,'exact values',flush=True)
reciprocal = median = t_checks = 0
for (d1,d2,p),x in values.items():
comp=(d2,d1,str(1-QQ(p)))
assert comp in values
residual=x*values[comp]-1
assert residual.is_finite() and residual.contains_zero()
assert residual.abs() < R(10)**(-97)
reciprocal += 1
if d1==d2 and p=='1/2':
assert bounds[(d1,d2,p)] == (1,1)
median+=1
if d1==1:
target=t.ppf((1+float(QQ(p)))/2,d2)**2
assert abs(target-float(x))/float(x)<5e-10
t_checks+=1
if d1==2:
theory=R(d2)/2*((1-R(QQ(p)))**(-R(2)/d2)-1)
assert x.overlaps(theory)
if d2==2:
z=R(QQ(p))**(R(2)/d1)
theory=2*z/(d1*(1-z))
assert x.overlaps(theory)
if d1==d2==1:
theory=(R.pi()*R(QQ(p))/2).tan()**2
assert x.overlaps(theory)
for d1 in range(1,11):
for d2 in range(1,11):
for a,b in zip(PROBABILITIES,PROBABILITIES[1:]):
assert bounds[(d1,d2,a)][1] < bounds[(d1,d2,b)][0]
# Equal-tail conversion at alpha=.05, using p=.025 and .975.
l,h=values[(d1,d2,'1/40')],values[(d1,d2,'39/40')]
mass=elementary_cdf(h,d1,d2)-elementary_cdf(l,d1,d2)
assert mass.overlaps(R(QQ(19)/20))
assert (mass-R(QQ(19)/20)).abs()<R(10)**(-97)
nist_count=0; nist_not_nearest=[]
for d2,line in enumerate(NIST_95.strip().splitlines(),1):
for d1,text in enumerate(line.split(),1):
x=values[(d1,d2,'19/20')]
difference=(x-R(QQ(Fraction(text)))).abs()
assert difference<R(QQ(1)/1000), (d1,d2,text,str(x))
if difference>R(QQ(1)/2000):
nist_not_nearest.append({'d1':d1,'d2':d2,'printed':text,'computed':str(x)})
nist_count+=1
# Exact agreement helper is appropriate for rationals, not ball == ball.
rational22={p:bounds[(2,2,p)][0] for p in PROBABILITIES}
assert not agrees_with(rational22,lambda p:QQ(p)/(1-QQ(p)))
summary={'entries':len(records),'exact':exact_count,'certified_decimals':len(records)-exact_count,
'independent_cdf_endpoint_checks':2*(len(records)-exact_count),'reciprocal_checks':reciprocal,
'equal_degree_medians':median,'student_t_squares_scipy':t_checks,'two_sided_conversions':100,
'nist_within_one_last_place':nist_count,'nist_not_correctly_rounded':nist_not_nearest,'scipy_quantile_checks':len(records),
'max_scipy_relative_difference':max_scipy,'measurement':measure({k:e['number'] for k,e in zip(bounds,records)})}
print(json.dumps(summary,indent=2),flush=True)
return summary
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
gen=FQuantiles(); start=time.monotonic()
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
if '--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1':
stored=_stored_entries(numberdb.table('T455'))
records=[{'params':{k:str(v) for k,v in p.items()},
'number':stored[','.join(str(p[k]) for k in gen.parameters)]}
for p in gen.enumerate()]
else:
records=[]; balls=[]; min_bits=None
for params in gen.enumerate():
value=gen.value(params,gen.digits); balls.append(value)
records.append({'params':{k:str(v) for k,v in params.items()},'number':to_text(value,gen.digits)})
if hasattr(value,'rad'):
interval=RealIntervalField(700)(value)
low,high=map(QQ,interval.endpoints())
ratio=(high-low)/low
quality=ratio.numerator().nbits()-ratio.denominator().nbits()+1
min_bits=min(min_bits,-quality) if min_bits is not None else -quality
assert not exactness(balls)
print('minimum retained relative bits:',min_bits,flush=True)
# JSON is emitted so the wrapper can preserve the exact written data.
print('ENTRIES_JSON '+json.dumps(records),flush=True)
summary=inspect_entries(records)
print('elapsed seconds:',round(time.monotonic()-start,2),flush=True)
return summary
if __name__=='__main__':
main()