Efficient Calculation of the Joint Distribution of Order Statistics

TL;DR

This paper introduces novel recursive methods for efficiently computing the joint distribution of order statistics of independent variables, enabling exact and rounded calculations, with applications in multiple hypothesis testing.

Contribution

It generalizes existing recursive formulas to improve numerical computation of joint distributions of order statistics, both exactly and approximately.

Findings

Developed generalized recursive formulas for joint distribution calculation.

Achieved exact and faithfully rounded numerical results.

Discussed applications in stepwise multiple hypothesis testing.

Abstract

We consider the problem of computing the joint distribution of order statistics of stochastically independent random variables in one- and two-group models. While recursive formulas for evaluating the joint cumulative distribution function of such order statistics exist in the literature for a longer time, their numerical implementation remains a challenging task. We tackle this task by presenting novel generalizations of known recursions which we utilize to obtain exact results (calculated in rational arithmetic) as well as faithfully rounded results. Finally, some applications in stepwise multiple hypothesis testing are discussed.

Tables

Table 1. Table 1 : Calculation of the probability that uniform order statistics are bounded above by b i ∈ [ 0 , 1 ] 1 ​ 1 subscript 𝑏 𝑖 superscript 0 1 1 1 b_{i}\in[0,1]^{1}1 where b i = 2 − 10 subscript 𝑏 𝑖 superscript 2 10 b_{i}=2^{-10} if i ≤ 10 𝑖 10 i\leq 10 and b 11 = 2 − 1 subscript 𝑏 11 superscript 2 1 b_{11}=2^{-1} . The rows give the intermediate steps of the algorithms. The first column reports the first few non-zero digits of the exact probabilities (computed in rational arithmetic), the second column reports the steps of Steck’s recursion (calculated in double precision floating point arithmetic) and the fourth column reports the steps of Bolshev’s recursion (also calculated in double precision floating point arithmetic). In the third and the fifth column the relative error of the intermediate value is reported.

	Exact Probability	Steck	Rel. Err. (Steck)	Bolshev	Rel. Err. (Bolshev)
2	9.76562E-04	9.76562E-04	0.00000E+00	9.76562E-04	0.00000E+00
3	9.53674E-07	9.53674E-07	0.00000E+00	9.53674E-07	0.00000E+00
4	9.31323E-10	9.31323E-10	0.00000E+00	9.31323E-10	0.00000E+00
5	9.09495E-13	9.09495E-13	0.00000E+00	9.09495E-13	0.00000E+00
6	8.88178E-16	8.88178E-16	0.00000E+00	8.88178E-16	0.00000E+00
7	8.67362E-19	8.67362E-19	0.00000E+00	1.73472E-18	1.00000E+00
8	8.47033E-22	8.47033E-22	0.00000E+00	1.10114E-20	1.20000E+01
9	8.27181E-25	8.27181E-25	0.00000E+00	6.85071E-21	8.28100E+03
10	8.07794E-28	8.07794E-28	0.00000E+00	-2.70517E-20	3.34884E+07
11	7.88861E-31	7.88861E-31	0.00000E+00	-1.12683E-16	1.42842E+14
12	4.33103E-30	-1.75898E-20	4.06134E+09	2.83880E-16	6.55456E+13

Table 2. Table 2 : Average Power of the Benjamini-Hochberg procedure (controlling the FDR at α = 0.05 𝛼 0.05 \alpha=0.05 ) for m 𝑚 m independent two-sided one sample z-tests (sample size N = 5 𝑁 5 N=5 , common variance σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 , μ 0 = 0 , μ 1 = 1 formulae-sequence subscript 𝜇 0 0 subscript 𝜇 1 1 \mu_{0}=0,\mu_{1}=1 , cf. ( 10 )) when m 0 subscript 𝑚 0 m_{0} hypotheses are true. The bold values are exactly those in (Glueck et al., 2008b , Table 2) .

$m$	$m_0=0$	$m_0=1$	$m_0=2$	$m_0=3$	$m_0=4$	$m_0=5$
2	0.56539	0.50342
3	0.54576	0.49842	0.44439
4	0.53446	0.49583	0.45256	0.40451
5	0.52712	0.49440	0.45819	0.41837	0.37494
6	0.52201	0.49357	0.46241	0.42840	0.39148	0.35175
7	0.51827	0.49310	0.46574	0.43606	0.40399	0.36955
8	0.51543	0.49285	0.46846	0.44214	0.41383	0.38349
9	0.51320	0.49273	0.47073	0.44710	0.42178	0.39473
10	0.51142	0.49270	0.47266	0.45124	0.42836	0.40398
11	0.50997	0.49272	0.47434	0.45475	0.43390	0.41174
12	0.50877	0.49278	0.47580	0.45777	0.43863	0.41833
13	0.50776	0.49286	0.47709	0.46039	0.44271	0.42401
14	0.50690	0.49295	0.47823	0.46268	0.44627	0.42895
15	0.50616	0.49306	0.47925	0.46472	0.44941	0.43329
16	0.50552	0.49316	0.48018	0.46653	0.45219	0.43712
17	0.50497	0.49327	0.48101	0.46815	0.45467	0.44053
18	0.50447	0.49338	0.48177	0.46962	0.45690	0.44358
19	0.50404	0.49348	0.48246	0.47094	0.45891	0.44634
20	0.50365	0.49358	0.48309	0.47215	0.46074	0.44883
21	0.50330	0.49368	0.48367	0.47325	0.46240	0.45109
22	0.50298	0.49378	0.48421	0.47427	0.46392	0.45316
23	0.50269	0.49387	0.48471	0.47520	0.46532	0.45505
24	0.50243	0.49395	0.48517	0.47605	0.46660	0.45679
25	0.50219	0.49404	0.48559	0.47685	0.46779	0.45840
26	0.50198	0.49412	0.48599	0.47759	0.46889	0.45988
27	0.50178	0.49420	0.48637	0.47828	0.46991	0.46126
28	0.50159	0.49427	0.48671	0.47892	0.47086	0.46254
29	0.50142	0.49434	0.48704	0.47952	0.47175	0.46373
30	0.50126	0.49441	0.48735	0.48008	0.47258	0.46485
31	0.50111	0.49447	0.48764	0.48060	0.47336	0.46589
32	0.50097	0.49453	0.48791	0.48110	0.47409	0.46687
33	0.50084	0.49459	0.48817	0.48157	0.47478	0.46779
34	0.50072	0.49465	0.48841	0.48201	0.47542	0.46866
35	0.50060	0.49470	0.48864	0.48242	0.47604	0.46948
36	0.50050	0.49475	0.48886	0.48282	0.47661	0.47025
37	0.50040	0.49480	0.48907	0.48319	0.47716	0.47098
38	0.50030	0.49485	0.48926	0.48354	0.47768	0.47167
39	0.50021	0.49489	0.48945	0.48388	0.47817	0.47233
40	0.50012	0.49494	0.48963	0.48420	0.47864	0.47295
41	0.50004	0.49498	0.48980	0.48451	0.47909	0.47354
42	0.49996	0.49502	0.48996	0.48480	0.47951	0.47411
43	0.49989	0.49506	0.49012	0.48508	0.47992	0.47465
44	0.49982	0.49509	0.49027	0.48534	0.48031	0.47516
45	0.49975	0.49513	0.49041	0.48559	0.48068	0.47565
46	0.49969	0.49516	0.49055	0.48584	0.48103	0.47612
47	0.49963	0.49520	0.49068	0.48607	0.48137	0.47657
48	0.49957	0.49523	0.49081	0.48629	0.48169	0.47700
49	0.49951	0.49526	0.49093	0.48651	0.48201	0.47741
50	0.49946	0.49529	0.49104	0.48672	0.48230	0.47781

Table 3. Table 3 : Faithfully rounded calculation of the λ 90 − limit-from subscript 𝜆 90 \lambda_{90}- power of the Benjamini Hochberg procedure (controlling the FDR at α = 0.15 𝛼 0.15 \alpha=0.15 ) when applied to m = 200 𝑚 200 m=200 test-statistics with chi-squared distributions F ν , 0 subscript 𝐹 𝜈 0 F_{\nu,0} under the null hypothesis and F ν , μ subscript 𝐹 𝜈 𝜇 F_{\nu,\mu} under the alternative where ν = 2 ​ n − 2 𝜈 2 𝑛 2 \nu=2n-2 , μ = n 2 ​ θ 𝜇 𝑛 2 𝜃 \mu=\sqrt{\frac{n}{2}}\theta . The results are under the RM-model where π 0 = 𝔼 ​ ( M m ) m subscript 𝜋 0 𝔼 subscript 𝑀 𝑚 𝑚 \pi_{0}=\frac{\mathbb{E}(M_{m})}{m} . For comparison the fourth column contains the Monte Carlo estimates (sample size 1000 1000 1000 ) given in (Izmirlian, 2018 , Table 3) . The following column gives the faithfully rounded power calculated using our method and the last column states the absolute difference between the previous two columns dived by the standard deviation of the Monte Carlo estimation (which was estimated using 300 replicates). As expected most Monte Carlo approximations are within one or two standard deviations of the faithfully rounded result. For the twelfth row no value is given since the estimated standard deviation was zero. This is not unexpected since the faithfully rounded result is ≈ 4.58 ⋅ 10 − 7 absent ⋅ 4.58 superscript 10 7 \approx 4.58\cdot 10^{-7} which is two orders magnitude smaller than ( 300 ⋅ 1000 ) − 1 superscript ⋅ 300 1000 1 (300\cdot 1000)^{-1} .

	Eff Sz. $\theta$	$𝔼\left(M_m\right)$	n	est. ${\lambda }_{90}$-pwr	${\lambda }_{90}$-pwr	Diff in std
1	0.60000	5	70	0.24900	0.26691	1.23987
2	0.60000	5	80	0.39600	0.39977	0.24081
3	0.60000	5	90	0.53800	0.53479	0.18527
4	0.60000	5	100	0.65700	0.65626	0.05057
5	0.60000	20	50	0.02800	0.02538	0.48379
6	0.60000	20	60	0.15700	0.13890	1.69096
7	0.60000	20	70	0.37800	0.36864	0.61103
8	0.60000	20	80	0.59900	0.62231	1.50514
9	0.60000	60	40	0.00200	0.00143	0.43439
10	0.60000	60	50	0.09900	0.08584	1.49443
11	0.60000	60	60	0.49200	0.49307	0.06522
12	0.60000	100	30	0.00000	0.00000	Inf
13	0.60000	100	40	0.00600	0.00658	0.22041
14	0.60000	100	50	0.27000	0.30726	2.80406
15	0.80000	5	40	0.25200	0.26037	0.59552
16	0.80000	5	50	0.50200	0.49870	0.19588
17	0.80000	5	60	0.70400	0.70951	0.39928
18	0.80000	20	30	0.03600	0.03969	0.60858
19	0.80000	20	40	0.35700	0.36732	0.64187
20	0.80000	60	20	0.00000	0.00004	0.18826
21	0.80000	60	30	0.14700	0.15775	0.88238
22	0.80000	100	20	0.00300	0.00013	8.59402
23	0.80000	100	30	0.50600	0.48563	1.42549
24	1.00000	5	30	0.39200	0.39421	0.14658
25	1.00000	20	20	0.04500	0.04534	0.05420
26	1.00000	20	30	0.61400	0.63660	1.40588
27	1.00000	60	20	0.22500	0.19941	2.00156
28	1.00000	100	20	0.58600	0.57569	0.64138

Equations

$$i_j$$

$$I_j$$

$${\Psi }_{n_1,n_2}^{G_1,G_2}(\mathbf{b}):=\mathbb{P}\left(X_{1:n}\le b_1,\dots ,X_{n:n}\le b_n\right),\mathbf{b}=(b_1,\dots ,b_n{)}^{\top }\in {\mathbb{R}}^n,$$

$$G_1(X_i)\sim \{\begin{array}{ll}\text{Uni}[0,1],& F_i=G_1,\\ G_2\circ G_1^{-1},& F_i=G_2,\end{array}$$

$$\Psi (i_1,i_2)$$

$$\Psi (m_1,m_2)$$

$$M_{k_1,k_2}^{(m_1,m_2)}:=\left(\genfrac{}{}{0ex}{}{m_1}{k_1}\right)\left(\genfrac{}{}{0ex}{}{m_2}{k_2}\right)(1-b_{k_1+k_2+1}{)}^{m_1-k_1}\cdot (1-F(b_{k_1+k_2+1}){)}^{m_2-k_2}.$$

$$M_{k_1,k_2}^{(m_1+1,m_2)}$$

$$M_{k_1,k_2}^{(m_1,m_2+1)}$$

$$\Psi (m_1,m_2)$$

$$M_{k_1,k_2}^{(m_1,m_2)}:=\left(\genfrac{}{}{0ex}{}{m_1}{k_1}\right)\left(\genfrac{}{}{0ex}{}{m_2}{k_2}\right){\left(b_{m_1+m_2}-b_{k_1+k_2+1}\right)}^{m_1-k_1}{\left(F(b_{m_1+m_2})-F(b_{k_1+k_2+1})\right)}^{m_2-k_2}.$$

$$a(k,j):=\left(\genfrac{}{}{0ex}{}{k}{j}\right)\text{ and }a(k,j)=\{\begin{array}{ll}1& k=j\\ \frac{j+1}{k-j}\times a(k,j+1)& j<k,\end{array}$$

$$M_{0,j}^{(m_1,m_2)}$$

$$M_{j,m_2}^{(m_1,m_2)}$$

$$M_{k_1+1,k_2-1}^{(m_1,m_2)}$$

$$Q_{i_1,i_2}(m)$$

$$M_{k_1,k_2}^{i_1,i_2}(m)$$

$$\Psi (i_1,i_2)$$

$$a^{(m),1}(j):=(b_m-b_{m-1}{)}^j\text{ and }{a}^{(m),2}(j):={\left(F(b_m)-F(b_{m-1})\right)}^j,$$

$$M_{k_1,k_2}^{i_1,i_2}(m)=\left(\genfrac{}{}{0ex}{}{i_1}{k_1}\right)\left(\genfrac{}{}{0ex}{}{i_2}{k_2}\right)\times a^{(m),1}(i_1-k_1)\times a^{(m),2}(i_2-k_2).$$

$$b_i$$

$$\text{\#Operations}=n+(n+1)+\sum _{k=2}^n\left[2+\sum _{j=1}^{k-1}6\right]=3n^2+n-1$$

$$x_i:=\{\begin{array}{ll}{b}_1& \text{if }i=1\\ b_i-b_{i-1}& \text{if }1<i\le n\\ F(b_1)& \text{if }i=n+1\\ F(b_i)-F(b_{i-1})& \text{if }n+1<i\le 2n\end{array}$$

$$\tilde{\Psi }(i_1,i_2)\in \Psi (i_1,i_2)\cdot \left((1-\epsilon {)}^{i_1+i_2},(1+\epsilon {)}^{i_1+i_2}\right),$$

$$p_i\sim \{\begin{array}{ll}\text{Uni}[0,1]& \text{if }1\le i\le m_0,\\ F& \text{if }{m}_0+1\le i\le m,\end{array}$$

$${\text{SU}}_{\mathbf{t}}(\mathbf{p})$$

$${\mathbb{P}}_{m,{\pi }_0,F}\left(V({\text{SU}}_{\mathbf{t}},\mathbf{p})=j,R({\text{SU}}_{\mathbf{t}},\mathbf{p})=k\right)$$

$$=\left(\genfrac{}{}{0ex}{}{m}{k}\right)\left(\genfrac{}{}{0ex}{}{k}{j}\right){\tilde{\pi }}_0^j(1-{\tilde{\pi }}_0{)}^{k-j}G(t_k{)}^k{\Psi }_{m-k,0}^{\text{Uni}[0,1],F}(1-G(t_m),\dots ,1-G(t_{k+1}))$$

$${\mathbb{P}}_{m,m_0,F}\left(V({\text{SU}}_{\mathbf{t}},\mathbf{p})=j,R({\text{SU}}_{\mathbf{t}},\mathbf{p})=k\right)$$

$$=\left(\genfrac{}{}{0ex}{}{m_0}{j}\right)\left(\genfrac{}{}{0ex}{}{m-m_0}{k-j}\right)t_k^j{\left(F(t_k)\right)}^{k-j}{\Psi }_{\begin{array}{c}m-k-(m_0-j),\\ m_0-j\end{array}}^{\text{Uni}[0,1],\bar{F}}(1-t_m,\dots ,1-t_{k+1}),$$

$$\mathrm{FDP}({\text{SU}}_{\mathbf{t}},\mathbf{p}):=\frac{V({\text{SU}}_{\mathbf{t}},\mathbf{p})}{R({\text{SU}}_{\mathbf{t}},\mathbf{p})\vee 1}.$$

$${\mathrm{Pow}}_{\text{avg}}({\text{SU}}_{\mathbf{t}}):=\mathbb{E}\left[\frac{R({\text{SU}}_{\mathbf{t}},\mathbf{p})-V({\text{SU}}_{\mathbf{t}},\mathbf{p})}{m-M_0}\right],$$

$${\mathrm{Pow}}_{\lambda }({\text{SU}}_{\mathbf{t}}):=\mathbb{P}\left(\frac{R({\text{SU}}_{\mathbf{t}},\mathbf{p})-V({\text{SU}}_{\mathbf{t}},\mathbf{p})}{m-M_0}\ge \lambda \right)$$

Full text

\floatsetup

[table]font=scriptsize,style=plaintop

1 Introduction

The joint distribution of order statistics $X_{1:n},\ldots,X_{n:n}$ of stochastically independent random variables $X_{1},\ldots,X_{n}$ plays a pivotal role in the theory of empirical processes and in nonparametric statistics; see, e. g., Shorack and Wellner (2009) and Dickhaus (2018). For instance, the exact finite-sample null distributions of classical goodness-of-fit tests like the Kolmogorov-Smirnov and the Cramér-von Mises test as well as those of modern ”higher criticism” goodness-of-fit tests rely on such joint distributions; cf. Gontscharuk et al. (2015), Gontscharuk et al. (2016a), and Finner and Gontscharuk (2018) for recent developments and further references. In simultaneous statistical inference, the joint distribution of ordered $p$-values is needed to analyze the type I and type II error behaviour of stepwise rejective multiple test procedures; cf. Chapter 5 of Dickhaus (2014).

In the case that $X_{1},\ldots,X_{n}$ are identically distributed (we refer to this case as a one-group model), classical recursive methods like Bolshev’s recursion, Noe’s recursion, and Steck’s recursion allow for computing the joint cumulative distribution function (cdf) of $X_{1:n},\ldots,X_{n:n}$ exactly; cf. Section 9.3 of Shorack and Wellner (2009). A generalization of Steck’s recursion to two-group models has been introduced by Blanchard et al. (2014). The other aforementioned recursions can be generalized in an analogous manner, as we will demonstrate in Section 3 of the present work.

While conceptually appealing, numerical properties of the aforementioned recursions are not well understood yet, and existing implementations into computer software often refer to rule-of-thumb-type upper bounds on $n$ such that the respective implementation is trustworthy. For example, Art B. Owen reports in his implementation of the two-sided version of Noe’s recursion in C (see https://www.stat.washington.edu/jaw/RESEARCH/SOFTWARE/BERKJONES/BJ-RBJ-C-Code/noe.c) that the recursion works well for $n\leq 1000$ but ”For larger n (eg 1800 or more) […] unexplained odd behavior.” Similarly, in the R Package mutoss (cf. Blanchard et al. (2010)) the following comment is made on the implementation of Bolshev’s recursion: ”Because of numerical issues $n$ should not be greater than 100.” Recently, Moscovich and Nadler (2017) introduced a computational method for one-group models. However, they do not consider the numerical accuracy of their approach rigorously.

In this work, we contribute to the analysis of the numerical accuracy and the computational complexity of existing approaches for computing the joint distribution of $X_{1:n},\ldots,X_{n:n}$ in a mathematically rigorous manner. Furthermore, we provide novel computational techniques for one- and two-group models which are guaranteed to provide accurate results for arbitrary sample size $n$. The rest of the material is structured as follows. In Section 2, we introduce the relevant quantities. The (generalized) recursions for one- and two-group models are provided in Section 3, together with a rigorous analysis of their computational complexities and their numerical properties. Our proposed exact computational methods rely on rational arithmetic (Section 4) and on faithful rounding (Section 5), respectively. Applications in multiple hypothesis testing are given in Section 6, and we conclude with a discussion in Section 7. Lengthy proofs as well as pseudo code for the considered algorithms are deferred to the Appendix.

2 Order Statistics

Throughout the following sections, we let $[n]:=\{1,2,\ldots,n\}$ for a natural number $n\in\mathbb{N}$. Consider stochastically independent, real-valued random variables $X_{1},\ldots,X_{n}$, which are all driven by the same probability measure $\mathbb{P}$. Let $I_{1}:=[n]$, and recursively define

[TABLE]

for $j\in[n]$. Then we call $X_{i_{1}},\ldots,X_{i_{n}}$ the order statistics of $X_{1},\ldots,X_{n}$, which we will denote by $X_{1:n},\ldots,X_{n:n}$ in the remainder. The random variable $X_{i:n}$ will be called the $i$-th order statistic of the random vector $(X_{1},\ldots,X_{n})^{\top}$. Let $F_{i}$ denote the marginal cdf of $X_{i}$ for $i\in[n]$. This paper will present methods for the quick and numerically stable calculation of

[TABLE]

assuming that $\forall i\in[n]:F_{i}\in\{G_{1},G_{2}\}$ where $G_{1},G_{2}$ are two continuous distribution functions on $\mathbb{R}$ and with $n_{i}=\left|\left\{j\in[n]:F_{j}=G_{i}\right\}\right|$ denoting the number of $X_{j}$’s distributed according to $G_{i}$, $i=1,2$. Since it holds that

[TABLE]

it follows that $\Psi_{n_{1},n_{2}}^{G_{1},G_{2}}=\Psi_{n_{1},n_{2}}^{\text{Uni}{[0,1]},F}\circ G_{1}$, where $F:=G_{2}\circ G_{1}^{-1}$. Therefore, it is sufficient to consider the calculation of $\Psi_{n_{1},n_{2}}^{\text{Uni}{[0,1]},F}(\mathbf{b})$ for an arbitrary continuous distribution function $F:[0,1]\rightarrow[0,1]$ and argument $\mathbf{b}\in[0,1]^{n}$. In the sequel, we suppress the dependence on $F$ and $\mathbf{b}$ notationally, and write $\Psi(n_{1},n_{2}):=\Psi_{n_{1},n_{2}}^{\text{Uni}{[0,1]},F}(\mathbf{b})$ for notational convenience.

As outlined in the Introduction, for $n_{2}=0$ there exist many well known recursions (see e. g. Section 9.3 of Shorack and Wellner (2009)) for computing $\Psi(n_{1},n_{2})$. There are also newer approaches based on numerical integration(see Moscovich et al. (2016)) or based on the Poisson process (see Moscovich and Nadler (2017)). Unfortunately, the former cannot be easily generalized to the case $0<n_{2}<n$, since the Lebesgue density of an order statistic is in general not piece-wise constant. The latter is very fast due to usage of the Fourier transform, but numerically unstable for small values of the $b_{i}$’s. This can for instance be demonstrated using the thresholds of the well-known linear step-up test (cf. Benjamini and Hochberg (1995)) for control of the false discovery rate (FDR); see Figure 1. Glueck et al. (2008a) proposed an algorithm with exponential $O(n^{n})$ complexity (cf. (Glueck et al., 2008a, Theorem 4.2)), resulting in a very high computational effort for moderate or large values of $n$. However, the method of Glueck et al. (2008a) can be used to compute $k$-variate marginal distributions for $k<<n$, because in such cases the complexity of their approach reduces to $O(n^{k})$.

Since we are mostly concerned with the full joint distribution, we extend the approach suggested by Blanchard et al. (2014) and provide generalizations of Bolshev’s and Noe’s recursions. We compare them to the generalization of Steck’s recursion proposed by Blanchard et al. (2014) and demonstrate that the Bolshev recursion is suitable for exact computations in rational arithmetic, whereas Noe’s recursion is numerically stable when computed in fixed-precision floating point arithmetic.

All our numerical calculations were performed on a Windows 7 machine with an Intel(R) Core(TM) i7-4790 CPU with 32 gigabytes of RAM.

3 The Generalized Recursions

Let $n:=n_{1}+n_{2}$, $n_{1},n_{2}\in\mathbb{N}$. Furthermore, let $X_{1},\ldots,X_{n_{1}}\sim\text{Uni}[0,1]$ and $X_{n_{1}+1},\ldots,X_{n_{1}+n_{2}}\sim F$ be jointly stochastically independent. Let for $0\leq i_{1}\leq n_{1}$ and $0\leq i_{2}\leq n_{2}$

[TABLE]

where $M:=[i_{1}]\bigcup\left\{n_{1}+j\left|j\in[i_{2}]\right.\right\}$, $[0]:=\emptyset$, $X_{i:M}$ denotes the $i$-th order statistic of $(X_{j})_{j\in M}$, and $(b_{i})_{i\in[n]}$ is an increasing sequence with values in $[0,1]$. The following subsections provide formulas for efficiently calculating $\Psi(n_{1},n_{2})$, and we discuss their computational and numerical properties.

3.1 Generalization of Bolshev’s Recursion

Lemma 1** (Generalization of Bolshev’s Recursion).**

The function $\Psi$ from (1) satisfies the recursion

[TABLE]

where

[TABLE]

Moreover, we have the following recursive relationships for $M$.

[TABLE]

For $n_{1}=0$ or $n_{2}=0$ this is simply the well-known Bolshev recursion.

3.2 Generalization of Steck’s Recursion

Lemma 2** (Generalization of Steck’s Recursion).**

Let $b_{0}:=0$. Then $\Psi$ from (1) satisfies the recursion

[TABLE]

where

[TABLE]

Letting

[TABLE]

we can write

[TABLE]

Furthermore, we have the following recursion for $M$.

[TABLE]

for $0\leq m_{1}\leq n_{1}$ and $0\leq m_{2}\leq n_{2}$.

Proof.

See (Blanchard et al., 2014, Proposition 1) ∎

3.3 Generalization of Noe’s Recursion

Lemma 3** (Generalization of Noe’s Recursion).**

Let $b_{0}:=0$, $Q_{0,0}(0):=1$, $Q_{i_{1},i_{2}}(1):=b_{1}^{i_{1}}\cdot F(b_{1})^{i_{2}}$ and for $m>1$

[TABLE]

for $0\leq i_{1}\leq n_{1},0\leq i_{2}\leq n_{2},m\leq i_{1}+i_{2}\leq n$.

Then the function $\Psi$ from (1) satisfies

[TABLE]

for $i_{1}\leq n_{1}$ and $i_{2}\leq n_{2}$.

Letting

[TABLE]

we can write

[TABLE]

3.4 Computational Complexity and Numerical Properties

The computational complexity (defined to be the number of elementary arithmetic operations on floating point numbers) of each of the aforementioned recursions is given by the following lemma.

Lemma 4**.**

The proposed recursions can be implemented using

Bolshev

$O(n_{1}^{2}n_{2}^{2})$**

Steck

$O\left(n_{1}^{2}n_{2}^{2}\log_{2}(n_{1}n_{2})\right)$**

Noe

$O\left(n_{1}^{2}n_{2}^{2}(n_{1}+n_{2})\right)$**

elementary arithmetic operations (addition, subtraction, multiplication, division) and $O(n_{1}n_{2})$ memory (assuming fixed-precision storage of all results).

The results of Lemma 4 suggest that Noe’s recursion might not be the best choice. However, for small values of the $b_{i}$’s, Bolshev’s recursion and Steck’s recursion are inherently numerically unstable. Consider for example $n_{1}=11,n_{2}=0$ and

[TABLE]

Then both recursions, when implemented in double precision floating point arithmetic, result in negative values and huge relative errors (cf. Table 1) which can be explained by inaccurate or catastrophic cancellation, respectively.

Noe’s recursion, if implemented in a reasonable manner, never results in negative values. Furthermore by (Jeannerod and Rump, 2018, Equation (3)) the relative error is bounded (if the coefficients are computed with a bounded relative error) since all summands are non-negative.

Remark 1**.**

Noe’s recursion can be easily parallelized since the $Q_{i_{1},i_{2}}(m)$ can be, for any fixed $m$, computed in parallel.

4 Exact Evaluation of Bolshev’s Recursion

If only elementary arithmetic operations are utilized and the number of such operations is not too large it is feasible to exactly evaluate expressions using rational arithmetic.111Our C++ implementation is based on The GNU Multiple Precision Arithmetic Library. We show that this is indeed the case for the Bolshev recursion as well its generalization for the two-group case presented in Section 3.1.

First we consider the case where $n_{2}=0$, hence $n_{1}=n$: Even though Bolshev’s recursion involves binomial coefficients our proposed Algorithm 1 (cf. the Appendix) for the one-group case evaluates it using only

[TABLE]

elementary arithmetic operations (addition, subtraction, multiplication, division). For an illustration, considering the sequence $b_{i}\equiv b_{i}^{(n)}:=0.05\times i/n$ we observed the execution times depicted in Figure 2.

Remark 2**.**

Our proposed Algorithm 2 (cf. the Appendix) implements the two-group case in $O(n_{1}^{2}n_{2}^{2})$ elementary arithmetic operations. Consequently, for equal sample sizes $n_{1}=n_{2}=\ell$ the number of operations is of $\operatorname{O}(\ell^{4})$. Notice that this is a marked improvement over the exponential complexity $\ell^{\ell}$ reported by (Glueck et al., 2008a, Theorem 4.2). Figure 3 illustrates the observed execution time for calculating $f(\ell):=\Psi(\ell,\ell)$ and $b_{i}\equiv b_{i}^{(n)}:=0.05\times i/n$, where $n=n_{1}+n_{2}=2\ell$.

For not necessarily equal sample sizes $n_{1}\geq 1$ and $n_{2}\geq 1$, our implementation of Algorithm 2 needs $f(n_{1},n_{2}):=1.5\cdot n_{1}^{2}\cdot n_{2}^{2}+4.5\cdot(n_{1}^{2}\cdot n_{2}+n_{1}\cdot n_{2}^{2})+3\cdot(n_{1}+n_{1}^{2}+n_{2}+n_{2}^{2})+7.5\cdot n_{1}\cdot n_{2}+2$ arithmetic operations.

Since the cdf $F$ of many interesting distributions is not available in a closed form the thresholds $F(b_{i})$ might either not be exactly calculable or simply not exactly representable as rational numbers. Lemma 5 analyzes the error propagation when $F$ and / or $\mathbf{b}$ are inexact.

Lemma 5**.**

Let

[TABLE]

and denote by $\left(\tilde{x}_{i}\right)_{i\in[n]}$ approximations thereof, which are obtained by replacing $\left(b_{i}\right)_{i\in[n]}$ and $F$ by approximations $(\tilde{b}_{i})_{i\in[n]}$ and $\tilde{F}$. If for $\varepsilon\in(0,1)$ it holds that for $\forall i\in[2n]:x_{i}\in(1-\varepsilon,1+\varepsilon)$ the it follows that for all $(i_{1},i_{2})\in[n_{1}]\times[n_{2}]$

[TABLE]

where $\tilde{\Psi}$ denotes the approximation of $\Psi$ obtained by using $\tilde{F}$ and $\tilde{b}_{i}$ instead of $F$ and $b_{i}$.

5 Faithfully Rounded Evaluation of Noe’s Recursion

We implemented faithfully rounded222That is, the result is either exact (if the exact value is a floating point number) or it is one of the two closest floating point numbers. floating-point computations as described by Rump and Lange (2017) as a portable single-header C++11 library.333Available at https://github.com/jvschroeder/PairArithmetic/. Utilizing this library we implemented the generalization of Noe’s recursion presented in Section 3.3 obtaining faithfully rounded results if no underflow occurs.444In our experience this is usually the case if the values of $\Psi$ are not too close to the smallest (in absolute value) normal double, which equals $2^{-1022}\approx 2.225\cdot 10^{-308}$ on most computer architectures. In case of an underflow the results are smaller than the true values of $\Psi$, but never less than zero. Parallelization was implemented using Intel®Threading Building Blocks (TBB).

Notice that Noe’s recursion (and our generalization thereof) satisfies the NIC principle (No Inaccurate Cancellation, cf. (Rump and Lange, 2017, Definition 2.2)), that is there are no sums where at least one summand is not an input to the algorithm and the summands have opposite signs. Thus, by examining the evaluation tree (cf. (Rump and Lange, 2017, Definition 2.2)) of a concrete implementation, it is possible to calculate a number $k=k(n_{1},n_{2})$ according to Equation (11) of Rump and Lange (2017). The result will be faithfully rounded if no under- or overflow occurs, and $k\leq 2^{26}-2$ (when utilising a double precision floating point numbers). For our concrete implementation we obtain $k(n_{1},n_{2})=n_{1}\cdot n_{2}+8\cdot(n_{1}+n_{2})-7$, provided that $n_{1}+n_{2}\geq 2$. Thus (assuming that no over- or underflow occurs) the result is guaranteed to be faithfully rounded if $n_{1},n_{2}\leq 8184$. Our implementation could be, in terms of $k$, significantly improved by using binary summation. For example, for $n_{1}=n_{2}=400$ we obtain $k(400,400)=166{,}398$, while the corresponding number of $k$ in the case of binary summation would equal $17{,}421$. The latter improvement however comes at an additional computational cost, and may be considered mostly of theoretical interest since the calculation for $n_{1}=n_{2}=400$ already takes approximately $27$ minutes on a $4$ core Intel CPU. Figure 3 compares the runtime of our implementation of our generalization of Noe’s recursion to that of the algorithm from the previous section. It becomes apparent that Noe’s recursion with faithful rounding is much faster then Bolshev’s recursion implemented in rational arithmetic. For practical applications, we therefore recommend Noe’s recursion with faithful rounding, at least if a fixed numerical precision is sufficient.

6 Applications in Multiple Hypothesis Testing

As discussed by Roquain and Villers (2011), the values of $\Psi(i_{1},i_{2})$ for all $0\leq i_{1}\leq n_{1}$ and $0\leq i_{2}\leq n_{2}$ are important building blocks for calculating the joint distribution of the number of rejections $R$and the number of false rejections $V$ for step-up multiple tests. The random variables $V$ and $R$ play an important role when analyzing the type I and type II error behavior of such multiple tests. One important observation is, that the previously discussed recursions for calculating $\Psi(n_{1},n_{2})$ also calculate all such $\Psi(i_{1},i_{2})$’s as intermediate results.

Following Blanchard et al. (2014) we consider $m\geq 2$ null hypotheses $H_{1},\ldots,H_{m}$ which are simultaneously under consideration under one and the same statistical model. We assume that associated $p$-values $p_{1},\ldots,p_{m}$ are available on which the multiple test operates. Furthermore, we assume that $p_{1},\ldots,p_{m}$ (regarded as random variables) are jointly distributed according to one of the following models

$\text{FM}(m,m_{0},F)$

The $p_{i}$’s are stochastically independent with marginal distributions

[TABLE]

where $m_{0}$ denotes the number of true null hypotheses among $H_{1},\ldots,H_{m}$ and $F$ is a given continuous cdf on $[0,1]$.

$\text{RM}(m,\pi_{0},F)$

Let $M_{0}$ denote a binomially distributed random variable, $M_{0}\sim\mathcal{B}(m,\pi_{0})$. Conditionally on $M_{0}=m_{0}$, the $p_{i}$’s are jointly distributed according to $\text{FM}(m,m_{0},F)$.

A multiple test operating on $\mathbf{p}=(p_{1},\ldots,p_{m})^{\top}$ is a measurable mapping $\varphi:[0,1]^{m}\rightarrow\mathcal{P}([m])$, where hypothesis $H_{i}$ is rejected iff $i\in\varphi(\mathbf{p})$. Under $\text{FM}(m,m_{0},F)$ denote by $M_{0}\equiv m_{0}$ a constant random variable. Then the (random) number of rejections of the multiple test $\varphi$ is given by $R(\varphi,\mathbf{p}):=\left|\varphi(\mathbf{p})\right|$, and $V(\varphi,\mathbf{p}):=\left|\varphi(\mathbf{p})\cap[M_{0}]\right|$ is the (random) number of false rejections (type I errors).

In the following we will consider step-up procedures $\varphi=\text{SU}_{\mathbf{t}}$ with critical values $\mathbf{t}=(t_{1},\ldots,t_{m})^{\top}\in(0,1)^{m}$ such that $t_{1}\leq\ldots\leq t_{m}$. The corresponding decision rule can be written as

[TABLE]

Summarizing results of Roquain and Villers (2011), the joint distribution of $R$ and $V$ for any step-up procedure $\text{SU}_{\mathbf{t}}$ has the following properties.

Lemma 6**.**

Let $0\leq j\leq k\leq m$.

1. (i)

Under the unconditional model $RM(m,\pi_{0},F)$ it holds that

[TABLE]

where $\tilde{\pi}_{0}:=\pi_{0}t_{k}/G(t_{k})$ and $G(t):=\pi_{0}t+(1-\pi_{0})F(t)$. 2. (ii)

Under the conditional model $FM(m,m_{0},F)$ it holds that

[TABLE]

where $0\leq j\leq m_{0}$ and $\bar{F}(t):=1-F(t)$.

Combining Lemma 6 with the previously discussed efficient evaluation of $\Psi$ it is possible to calculate various summary statistics pertaining to the joint distribution of $(V,R,M_{0})$ under the above models.

Definition 1**.**

**

• (a)

The FDR of $\text{SU}_{\mathbf{t}}(\mathbf{p})$ is given by the expectation of the false discovery proportion (FDP) of $\text{SU}_{\mathbf{t}}(\mathbf{p})$, which is given by

[TABLE]

• (b)

Considering the number of correct rejections $R(\text{SU}_{\mathbf{t}},\mathbf{p})-V(\text{SU}_{\mathbf{t}},\mathbf{p})$ the average power of $\text{SU}_{\mathbf{t}}(\mathbf{p})$ is given by

[TABLE]

where the convention $\frac{0}{0}=0$ is utilized and where $\mathbb{E}=\mathbb{E}_{m,\pi_{0},F}$ (under $RM(m,\pi_{0},F)$) or $\mathbb{E}=\mathbb{E}_{m,m_{0},F}$ (under $FM(m,m_{0},F)$), respectively.

• (c)

*The *$\lambda-$power is the probability of rejecting at least $\lambda\cdot(m-M_{0})$ of the false hypotheses:

[TABLE]

where, again, the convention $\frac{0}{0}=0$ is utilised and where $\mathbb{P}=\mathbb{P}_{m,\pi_{0},F}$ (under $RM(m,\pi_{0},F)$) or $\mathbb{P}=\mathbb{P}_{m,m_{0},F}$ (under $FM(m,m_{0},F)$), respectively.

In order to provide some numerical illustrations, we first consider the average power (cf. (8)) under $FM(m,m_{0},F)$ where

[TABLE]

and $N=5$. This is the setting considered in (Glueck et al., 2008b, Table 2) where the average power for $m\leq 5$ was calculated for the Benjamini-Hochberg procedure (controlling the FDR at $\alpha=0.05$) for $m$ independent two-sided one sample z-tests. In our notation, the Benjamini-Hochberg (linear step-up test) procedure equals $\text{SU}_{\mathbf{t}}$ with $t_{i}=i\alpha/m$ for $i\in[m]$. Table 2 illustrates the results obtained for $m,m_{0}\leq 50$. Due to space constraints only the first six columns (corresponding to $m_{0}\leq 5$) are presented. The calculation of the full table (not presented here) took less than a second for an $m$ one magnitude larger (50 instead of 5) than the one considered by Glueck et al. (2008b). Figure 4 illustrates the time needed to calculate one row of such a table corresponding to some $m\in\mathbb{N}$ when utilizing our proposed algorithms.

As a second example, we consider the computation of $\operatorname{Pow}_{\lambda}(\text{SU}_{\mathbf{t}})$ from (9). Again, we choose $\mathbf{t}$ as in the Benjamini-Hochberg case. An asymptotic approximation of this quantity for

[TABLE]

where $F_{\nu,\mu}$ denotes the distribution function of a non-central chi-squared random variable with $\nu$ degrees of freedom and non-centrality parameter $\mu$, is given in (Izmirlian, 2018, Table 3). Our results can be used to calculate $\operatorname{Pow}_{\lambda}(\text{SU}_{\mathbf{t}})$. Table 3 gives the faithfully rounded values for the $0.9$-power for the parameters considered in (Izmirlian, 2018, Table 3).

We conclude by giving an example for the exact distribution of the FDP which shows why the FDR is not always an appropriate summary statistic. Consider again the multiple two-sided z-test described in Glueck et al. (2008b), that is $F$ given by (10), for $N=m=50$ and $m_{0}=5$. It is clear that in Figure 5 the distribution of the FDP is neither symmetric about its mean (the FDR, which is depicted as dotted vertical line) nor concentrated around the FDR. A similar argumentation has been used by, among others, Blanchard et al. (2014) and Delattre and Roquain (2015) in order to motivate the computation of the full distribution of the FDP and to control its quantiles. The latter task is inherently computationally demanding.

7 Discussion

We have presented computationally efficient and numerically stable methods for calculating the joint distribution of order statistics. Such joint distributions have a multitude of important applications that require their repeated evaluation (to numerically solve optimization problems). Apart from the applications that we have presented in Section 6, they include, among others, the calibration goodness-of-fit tests with equal local levels (see Section 1.4 of Gontscharuk et al. (2016b)) and the adjustment of the asymptotically optimal rejection curve as proposed by Finner et al. (2012), see Equation (19) in their paper, and (Finner et al., 2009, Equation (6.1)) with the goal of obtaining valid critical values for a step-up-down procedure (guaranteeing strict FDR control). The latter applications have not been considered explicitly in the present work, because they merely refer to the one-group case. For this case, the methods of Moscovich et al. (2016) are already sufficiently accurate and fast.

Future extensions to our methods could include a normalization of the exponents (in Noe’s recursion) to avoid underflows and the exploration of potential efficiency gains in the exact computation of Bolshev’s recursion by a trade-off between the memory consumption and the frequency of normalizations of the intermediate rational numbers.

A preliminary version of our planned package (which utilizes RCPP, cf. Eddelbuettel (2013)) for the R language (R Core Team (2017)) is available at https://github.com/jvschroeder/OrdStat/ and can be installed using the devtools package:

Ψinstall.pacakges("devtools") Ψdevtools::install_github("jvschroeder/OrdStat")

The code used to generate the graphics and numerical examples is available at https://github.com/jvschroeder/OrdStatExamples/. The graphics were created using ggplot2 (Wickham (2016)) and the R package tikzDevice.

Acknowledgments

Jonathan von Schroeder is supported by the Deutsche Forschungsgemeinschaft (DFG) within the framework of RTG 2224 ”$\pi^{3}$: Parameter Identification - Analysis, Algorithms, Applications”.

Appendix A Proofs

• Proof of Lemma 1.

Let $c_{X}(t):=\left|\left\{U_{i}\leq t|i\in X\right\}\right|$, $m_{1}\in[n_{1}]$, $m_{2}\in[n_{2}]$ and $m:=m_{1}+m_{2}$. Then, mimicking the approach of (Shorack and Wellner, 2009, p. 367 ff.) and (Blanchard et al., 2014, Proposition 1), it holds that

[TABLE]

Since this holds for any $m_{1}\in[n_{1}]$, $m_{2}\in[n_{2}]$ it follows that

[TABLE]

where $M$ is given by (2).

The recursions for $M$ follow from the definition of the binomial coefficient and routine calculations. ∎

• Proof of Lemma 3.

Let $(a_{i})_{i\in[n]}$ be an increasing sequence in $[0,1]$ such that $\forall i\in[n]:a_{i}<b_{i}$. Using notation similar to that of (Shorack and Wellner, 2009, p. 362 ff.) let $i:=i_{1}+i_{2}$ and

[TABLE]

where the $c_{j}$ are the $2i$ boundaries $a_{1},\cdots,a_{i},b_{1},\cdots,b_{i}$ arranged in any ascending order. To extend the induction in (Shorack and Wellner, 2009, p. 364 ff.) to the two-group case we only need to mimick the the approach of (Blanchard et al., 2014, Proposition 1) to provide a recursive formula for $Q_{i_{1},i_{2}}(m)$, given that $h(m+1)-1\leq i\leq g(m-1)$ (where $g,h$ are given by (Shorack and Wellner, 2009, p. 362, Eq. (17) and Eq. (18))). To this end note that in this case

[TABLE]

and $\bar{M}:=[i]\setminus M$ denotes the complement of $M$.

If all $a_{i}=0$, then it holds that

[TABLE]

which implies

[TABLE]

From the (11) it follows that $Q_{i}(n)=\begin{cases}1&i=0\\ 0&i>0\end{cases}$. Thus ${Q_{i_{1},i_{2}}\left(n+1\right)=b_{1}^{i_{1}}\cdot F(b_{1})^{i_{2}}}$ and for $m>n+1$ it holds that

[TABLE]

which needs to be calculated for $m-n\leq i\leq n$.

Let $b_{0}:=0$, $\tilde{Q}_{0,0}(0):=1$, $\tilde{Q}_{i_{1},i_{2}}(1):=b_{1}^{i_{1}}\cdot F(b_{1})^{i_{2}}$ and for $m>1$

[TABLE]

(which needs to be calculated for $\forall(i_{1},i_{2})\in[n_{1}]\times[n_{2}]:m\leq i_{1}+i_{2}\leq n$). Then it holds that

[TABLE]

for all $i_{1}\in[n_{1}],i_{2}\in[n_{2}]$. ∎

• Proof of Lemma 4.

Counting the numer of operations in the loops of Algorithm 2 it follows that

[TABLE]

holds. For the space complexity simply note that, to use the recursions for $M^{(m_{1},m_{2})}_{k_{1},k_{2}}$, we need to keep track of at most $M^{(m_{1},m_{2}-1)}_{k_{1},k_{2}}$ and $M^{(m_{1}-1,m_{2})}_{k_{1},k_{2}}$ (which is pessimistic - cf. algorithm 2).

For Steck’s recursion first note that, using exponentiation by squaring, one can calculate $a^{n}$ in $O(\log_{2}(n))$ multiplications (where $n\in\mathbb{N}$, c.f. (Knuth, 1981, p. 442, Algorithm A)). Thus the first row and last column of $M^{(m_{1},m_{2})}$ (cf. equations (4) (5)) can be calculated in

[TABLE]

Thus to calculate all the coefficient matrices we need at most

[TABLE]

arithmetic operations (since we can calculate $a(m_{1},j_{1}),a(m_{2},j_{2})$ for $j_{1}\in[m_{1}],j_{2}\in[m_{2}]$ in $O(m_{1}+m_{2})$ using (2)). It remains to note that (3) needs at most $O(m_{1}m_{2})$ arithmetic operations. For the space complexity simply note that we do not need to keep track of the previous coefficient matrices.

For Noe’s recursion first note that, for every $m\in\mathbb{N}$, we can calculate $a^{(m),1}(j_{1}),a^{(m),2}(j_{2})$ for $j_{1}\in[i_{1}],j_{2}\in[i_{2}]$ in $O(j_{1}+j_{2})$. Furthermore, using (2), the binomial coefficients in (6) can be calculated in $O(i_{1}+i_{2})$. Thus $M^{i_{1},i_{2}}$ can be calculated in $O(i_{1}i_{2})$. Thus $Q_{i_{1},i_{2}}(m)$ (assuming the necessary $Q_{\cdot,\cdot}(m-1)$ have already been calculated) is $O(i_{1}i_{2})$. Therefore the overall computational complexity is at most

[TABLE]

For the space complexity simply note, again, that we do not need to keep track of the previous coefficient matrices. ∎

• Proof of Lemma 5.

By Noe’s recursion (cf. lemma 3) the probability $\Psi(i_{1},i_{2})$ can be obtained by evaluating a polynomial of degree $i_{1}+i_{2}$ with only positive coefficients at $\mathbf{x}\in\mathbb{R}^{2n}$ (where $\mathbf{x}$ is given by (7)). It is therefore sufficient to show that the statement is true for such polynomials when applied to non-negative arguments. To provide a concise proof we utilise interval arithmetic (cf. Kearfott (1996)). Due to linearity (since all coefficients and inputs are non-negative) it is sufficient to show the claim for monomials $p(x_{1},\cdots,x_{2n}):=\prod_{i=1}^{2n}x_{i}^{a_{i}}$ with $a_{i}\in\mathbb{N}$, $\sum_{i=1}^{2n}a_{i}\leq i_{1}+i_{2}$. Since $0<1-\varepsilon<1$ and $1+\varepsilon>1$ it follows that

[TABLE]

which implies

[TABLE]

due to (Kearfott, 1996, Equation (4)). ∎

• Proof of Lemma 6.

By (Roquain and Villers, 2011, Theorem 3.1) under the unconditional model $RM(m,\pi_{0},F)$ and for a step-up procedure $\text{SU}_{\mathbf{t}}$:

[TABLE]

where

[TABLE]

and $\mathbb{P}$ denotes $\mathbb{P}_{m,\pi_{0},F}$.

Furthermore under the conditional model $FM(m,m_{0},F)$ and for a step-up procedure $\text{SU}_{\mathbf{t}}$ it holds (by (Roquain and Villers, 2011, Section 5.3)) that (where $\bar{F}(t):=1-F(1-t)$):

[TABLE]

where $\mathbb{P}$ denotes $\mathbb{P}_{m,m_{0},F}$. ∎

Appendix B Algorithms