back to table · edit · history · where entries came from · files · download
6516 bytes, as of the version from 2026-10-28 22:35 (current). Recorded here, not run.
# Studentized range quantiles: definitions and certification
This attachment belongs to numberdb.org/T457.
Let Z_1,...,Z_k be independent N(0,1) variables and V an independent
chi-squared variable with nu positive integer degrees of freedom. Write
Q=(max Z_i-min Z_i)/sqrt(V/nu), and store its lower-tail p-quantile.
The parameter order is nu,k,p. The computed rectangle has k=2,3,
nu=1,...,10,20,30,60,120, and the fourteen probabilities in generate.py.
The two-mean case is the ordinary absolute t comparison; three means is
the first nontrivial simultaneous comparison. The limited number of means
makes 100-digit, certified values possible with the method described here.
This is a partial table of the full family, not a restriction of its domain.
## Reduction for two means
(Z_1-Z_2)/sqrt(2) is standard normal, independently of V. Hence
Q_{2,nu}=sqrt(2)|T_nu| in distribution, and
q_{2,nu,p}=sqrt(2) t_{nu,(1+p)/2}.
The CDF at q>0 is I_x(1/2,nu/2), x=q^2/(2nu+q^2).
The generator evaluates I with Arb's hypergeometric function and the
complement relation, using the argument at most 1/2.
## Reduction for three means
Orthogonally project the Gaussian vector (Z_1,Z_2,Z_3) onto the plane
perpendicular to (1,1,1). Its two coordinates are independent standard
normals. The polar radius R has R^2 distributed as chi-squared with two
degrees of freedom, and its angle is uniform and independent of R.
The six oriented difference vectors e_i-e_j have norm sqrt(2) and equally
spaced directions in that plane. The range is their largest scalar product
with the residual vector. Within each sector of angle pi/3 it is
sqrt(2) R cos(theta), with -pi/6<=theta<=pi/6. The range is unchanged by
subtracting the sample mean. Conditioning on theta and V therefore gives
F_{3,nu}(q)=1-(6/pi) integral_0^(pi/6)
(1+q^2/(2nu cos(theta)^2))^(-nu/2) dtheta.
Here E(exp(-s V))=(1+2s)^(-nu/2), the chi-squared Laplace transform.
Every integrand is strictly decreasing in q>0, so this CDF is strictly
increasing from 0 to 1. Thus a strict sign change around p certifies the
unique quantile.
## A certified convergent series
Put a=q^2/(2nu), c=1/(1+a), b=(nu-1)/2 and d=nu/2. Substituting
y=sin(theta), with 0<=y<=1/2, rewrites the integral as
(1+a)^(-d) integral_0^(1/2) (1-y^2)^b (1-c y^2)^(-d) dy.
Let (1-z)^b (1-cz)^(-d)=sum A_n z^n. Logarithmic differentiation gives
A_0=1, A_{-1}=0,
(n+1)A_{n+1}=((1+c)n-b+dc)A_n+c(b-d-n+1)A_{n-1}.
The implementation recurs on A_n/4^n to keep intermediate values moderate.
It integrates N terms, obtaining (1/2) sum_{n=0}^{N-1} A_n/(4^n(2n+1)).
On the complex circle |z|=1/2 the generating function is analytic and
its absolute value is at most (3/2)^b 2^d <=3^nu. Cauchy's coefficient
bound gives |A_n|<=3^nu 2^n. Consequently the absolute integral tail is
at most 3^nu 2^(-N)/(2N+1), and in particular less than 3^nu 2^(-N).
The implementation adds the latter, slightly larger, error before scaling
by (1+a)^(-d) and 6/pi. It uses N=bits+2nu+32, so the tail is smaller than
2^(-bits-32). All coefficient operations and transcendental scaling use Arb
ball arithmetic, so rounding errors are included separately from this tail.
## Root certification and precision
mpmath quadrature and a secant iteration at 135 decimal digits propose a
candidate only. The proposal is converted to exact rational endpoints with
width 4*10^(-115). Arb at numberdb.bits(100, losing=192) bits, with explicit
escalation if necessary, must certify F(low)<p<F(high). A ball containing
that entire rational bracket is returned. Thus neither a midpoint rounded
into an interval nor a quadrature's estimated error is claimed as a proof.
For independent checks, checks.py uses elementary integration recurrences
for the Student density and for the three-mean angular integral. These
recurrences share neither the beta evaluation nor the infinite coefficient
recurrence used by the generator. They check opposite CDF signs at the
endpoints represented by every stored decimal, not merely its midpoint.
Independent controls precede those tests. SciPy's general two-dimensional
CDF provides a lower-precision check of the reduction to the angular
integral, rather than certifying any of the final 100 digits.
## Independent finite integration formulas
For two means put theta=atan(q/sqrt(2nu)). Then
F_{2,nu}(q)=C_{nu-1}(theta)/C_{nu-1}(pi/2), where
C_0(theta)=theta, C_1(theta)=sin(theta), and
C_j(theta)=sin(theta)cos(theta)^(j-1)/j+(j-1)C_{j-2}(theta)/j.
For three means write A=1+a, a=q^2/(2nu), and call its angular integral I.
If nu=2m+1, put u=asin(1/(2sqrt(A))). Substituting sin(theta)=sqrt(A)sin(u)
gives
I=A^(-m) sum_{j=0}^m binomial(m,j)(-a)^j T_j(u),
where T_0=u and T_j=tan(u)^(2j-1)/(2j-1)-T_{j-1}.
If nu=2m, substituting t=tan(theta) gives
I_m=integral_0^(1/sqrt(3)) dt/((1+t^2)(A+a t^2)^m).
The exact partial-fraction identity
1/((1+t^2)(A+a t^2)^m)
=1/((1+t^2)(A+a t^2)^(m-1))-a/(A+a t^2)^m
gives I_m=I_{m-1}-a K_m, with I_0=pi/6. For
u=atan(sqrt(a/(3A))), K_j=sqrt(A/a) A^(-j) C_{2j-2}(u).
Thus the CDF is directly the positive finite sum
F_{3,2m}(q)=(6/pi)sqrt(aA) sum_{j=1}^m A^(-j) C_{2j-2}(u).
These formulas use only a finite number of exact arithmetic and elementary
ball operations. The checker uses 2048 bits to leave ample room for cancellation in the
odd-nu tangent-power sum; it tests strict signs and finite enclosures
rather than accepting overlap alone. Three closed-form
controls return 1/2, and a fourth compares the angular formula with SciPy's
general CDF. Every row also rejects an interval shifted beyond its upper
endpoint. The one-degree-of-freedom closed forms and the two-mean,
two-degree-of-freedom closed form are checked on all relevant rows.
The 56 values in reference-values.json are the selected two- and three-mean
rows from the University of Washington's printed upper 0.05 and 0.01 table:
https://faculty.washington.edu/heagerty/Books/Biostatistics/TABLES/StudentizedRange/index.html
They are checks at the source's displayed precision, not a source of the
100-digit values. Fifty-five values agree within half a last-place unit. The remaining value,
nu=1,k=2,p=99/100, is printed as 90.03, whereas the exact Cauchy formula
sqrt(2)cot(pi/200) gives 90.02422668900665..., which rounds to 90.02.
The discrepancy is 0.5773311 printed units and is explicitly recorded, not
silently accepted as correct rounding. All 56 agree within one last-place
unit. The exact formula and independent ball CDF settle which value is right.