A Numerical Approach to Sequential Multi-Hypothesis Testing for Bernoulli Model

The problem of testing multiple hypotheses is one of the oldest problems in the sequential analysis.

A traditional approach to this problem is Bayesian. It is based on the assumption that the hypotheses come up with some probabilities called a priori (see Blackwell and Girshick 1954; Baum and Veeravalli 1994; Tartakovsky, Nikiforov, and Basseville 2015, among many others).

Despite that the optimal Bayesian solution can be characterized on the basis of general principles like dynamic programming or the theory of optimal stopping (Shiryaev 1978; Chow, Robbins, and Siegmund 1971), at least theoretically, there seems to exist a strong belief that the theoretical solution is too complex to be useful for practical purposes (see, for example Baum and Veeravalli 1994; Tartakovsky, Nikiforov, and Basseville 2015). An exception is the case of two simple hypotheses where the solution is given by the classical sequential probability ratio test (Wald’s SPRT, see Wald and Wolfowitz 1948).

For these reasons, approximate solutions of the problem have been proposed. One of the widely used tests, due to its simplicity, is the matrix sequential probability ratio test (MSPRT) by Armitage (1950). Tartakovsky, Nikiforov, and Basseville (2015) showed that the MSPRT is asymptotically optimal, as error probabilities go to 0.

Another approach that has received considerable attention through the decades is the so-called Kiefer-Weiss problem, consisting in the minimization of the maximum value of the expected sample number (ESS), over all possible parameter points (Kiefer and Weiss 1957). Lorden (1980) showed that the Kiefer-Weiss problem can be reduced to the minimization of the ESS evaluated at a specific parameter point, different from the hypothesized parameter values (so-called modified Kiefer-Weiss problem), and (in essence) used the method of Lagrange multipliers to characterize the solutions to the modified Kiefer-Weiss problem.

A generalization of the Kiefer-Weiss problem to the case of multiple hypotheses has been formulated in Tartakovsky, Nikiforov, and Basseville (2015) (Section 5.3) and received an asymptotic treatment in Section 5.3.1, ibid.

In this paper, we propose an approach to the optimal multi-hypothesis testing based on minimization of the weighted ESS evaluated at parameter points not necessarily coinciding with the hypothesized values, and then use the method of the Lagrange multipliers to reduce to the minimization of the Lagrangian function. Depending on the choice of the points for evaluating the ESS in the Lagrangian function, we obtain, in particular, the Bayesian and the Kiefer-Weiss settings, and more.

We apply the method of Novikov (2009b) and characterize the sequential tests minimizing the Lagrangian function, for any choice of multipliers. For practical applications, we propose the use of numerical methods for the Lagrange minimization, the evaluation of the characteristics (the error probabilities, the ESS, etc.), and for finding the multiplier values to comply with the restrictions on the error probabilities.

We illustrate the proposed methods in the particular case of sampling from a Bernoulli population, where we develop a complete set of computer algorithms for all the numerical tasks described above and implement them in the R programming language. The program code is available in a public GitHub repository in Novikov (2023).

Using the developed software, we run a series of numerical comparisons related to optimal properties of sequential multi-hypothesis tests in the Bernoulli model.

First, we evaluate the performance characteristics error of the MSPRT for a particular case of three hypotheses. The MSPRT is known to be asymptotically optimal, as the error probabilities go to 0, so the evaluations we carry out give an idea of how small the error probabilities should be in order that the asymptotic formulas for the ESS give a reasonably good approximation to the calculated values. We use $N=4000$ which, apparently, is sufficient for good approximations of the characteristics of non-truncated MSPRTs.

Other comparison we carry out is also related with the MSPRT. For a number of error probability levels, we numerically find both MSPRT and the optimal Bayes test (for uniform a priori distribution) matching the given error probabilities (up to some precision). The results show a very high efficiency of the MSPRT.

Also we propose a method for solving a multi-hypothesis version of the Kiefer-Weiss problem, and give a numerical example.

In Section 2, we adapt the results of Novikov (2009b) to the problem of minimization of weighted ESS calculated at arbitrary parameter points. In Section 3, we derive computational formulas for the Bernoulli model. Numerical results are presented in Section 4. Section 5 is a brief list of the results and suggestions for further work.

In this section, we formulate some settings for the problem of optimal multi-hypothesis testing and use the general results of Novikov (2009b) for characterisation of the respective optimal solutions.

We assume that independent and identically distributed (i.i.d.) observations $X_{1},X_{2},\dots,X_{n},\ \dots$ are potentially available to the statistician on the one-by-one basis, providing us with information about the unknown distribution of the data. Let us denote it $P_{\theta}$, where $\theta$ is some parameter identifying the distribution in a unique manner. We are concerned with the problem of distinguishing between a finite number of simple hypotheses $H_{1}:\theta=\theta_{1}$, $H_{2}:\theta=\theta_{2}$, $\dots$, $H_{k}:\theta=\theta_{k}$, $k\geq 2$.

We follow Novikov (2009b) in the notation and general assumptions.

In particular, we consider sequential multi-hypothesis test as a pair $\langle\psi,\phi\rangle$ of a stopping rule $\psi=(\psi_{1},\psi_{2},\dots,\psi_{n},\dots)$, and a (terminal) decision rule $\phi=(\phi_{1},\phi_{2},\dots,\phi_{n},\dots)$.

The elements of the stopping rule $\psi_{n}=\psi_{n}(x_{1},\dots,x_{n})$ are measurable functions taking values in $[0,1]$, where the value at $(x_{1},\dots,x_{n})$ is interpreted as the conditional probability, given the observations, to stop (randomization at the stopping time).

The elements of the decision rule $\phi_{n}=\phi_{n}(x_{1},\dots,x_{n})$ are measurable functions of observations such that $\phi_{n}=(\phi_{n}^{1},\dots,\phi_{n}^{k})$, and $\phi_{n}^{j}\geq 0$ and $\sum_{j=1}^{k}\phi_{n}^{j}(x_{1},\dots,x_{n})\equiv 1$. Given the data $(x_{1},\dots,x_{n})$ observed, $\phi_{n}^{j}(x_{1},\dots,x_{n})$ is interpreted as a conditional probability to accept hypothesis $H_{j}$, $j=1,\dots,k$ (randomization at the decision time).

The sequential test starts with observing $X_{1}=x_{1}$ (stage $n=1$). At each stage $n=1,2,\dots$ the test procedure stops with probability $\psi_{n}(x_{1},\dots,x_{n})$, given that $X_{1}=x_{1},\dots,X_{n}=x_{n}$ are observed, and proceeds to taking a terminal decision. If it does not stop, the test proceeds to taking one additional observation $X_{n+1}=x_{n+1}$ and going to stage $n+1$, etc., until the process eventually stops. When the test stops at any stage $n$ (this $n$ is called stopping time), a terminal decision is taken accepting hypothesis $H_{j}$ with probability $\phi_{n}^{j}(x_{1},\dots,x_{n})$, conditionally on $(x_{1},\dots,x_{n})$. Let us denote $\tau_{\psi}$ the stopping time (as a random variable) generated by the described process.

Let

($s_{1}^{\psi}(x_{1})=\psi_{1}(x_{1})$ by definition).

Then the expected sample size (ESS) of the test procedure is defined as

provided that $\sum_{n=1}^{\infty}E_{\theta}s_{n}^{\psi}=1$, - otherwize it is infinite by definition. Here and throughout the paper, $E_{\theta}$ is the symbol of mathematical expectation with respect to $P_{\theta}$. Also we use $s_{n}^{\psi}$ (without arguments) both for $s_{n}^{\psi}(x_{1},x_{2},\dots,x_{n})$ and for $s_{n}^{\psi}(X_{1},X_{2},\dots,X_{n})$, depending on the context. So do we when dealing with other functions like $\psi_{n}$, $\phi_{n}$, etc.

Other characteristics of a sequential test $\langle\psi,\phi\rangle$ are the error probabilities defined as

Another natural way to define error probabilities is less detailed:

In the case of two hypotheses the definitions are equivalent.

For $k=2$, the classical result of Wald and Wolfowitz (1948) states that the sequential probability ratio test (SPRT) minimizes both $E_{\theta_{1}}\tau_{\psi}$ and $E_{\theta_{2}}\tau_{\psi}$ in the class of sequential tests $\langle\psi,\phi\rangle$ such that

where $\alpha_{1}$ and $\alpha_{2}$ are the error probabilities of the SPRT.

To the best of our knowledge, no direct generalizations of this result exist for $k>2$. For this reason, we propose weaker settings.

Let us choose some parameter points $\vartheta_{i}$, $i=1,\dots,K$ and the weights $\gamma_{i}$, $i=1,\dots,K$ being these non-negative numbers such that $\sum_{i=1}^{K}\gamma_{i}=1$, $K\geq 1$. Formally, we propose to minimize the weighted ESS

over all sequential multi-hypothesis tests subject to

or to

where $\alpha_{ij}$ and $\alpha_{i}$ are some positive numbers.

To support this formulation, let us refer to a very practical context of optimal group-sequential testing in the case of two hypotheses. For testing the mean of a normal distribution with known variance, Eales and Jennison (1992) considered five settings for the ESS minimization under restrictions on the error probabilities. Four of them, namely, $F_{1}$ to $F_{4}$ (see Eales and Jennison 1992) are of type (1), with different choices of $K$, $\vartheta_{i}$ and $\gamma_{i}$. $F_{5}$ is also a kind of weighted ESS but of continuous type, which is quite possible to be treated by our method, but for the time being stays beyond the scope. Generalizations of these settings to the case of more than two hypotheses and infinite horizons are straightforward.

Given that the formulated problem is a problem of a minimization under restrictions, we want to use the Lagrange multipliers method. By the principle of the Lagrange method, to minimize $C_{\gamma,\vartheta}$ under restrictions (2) one should be able to minimize the Lagrangian function

with any constant multipliers $\lambda_{ij}\geq 0$, and to find the values of the multipliers for which equalities in (2) hold. Respectively, the problem of minimization under conditions (3) reduces to minimization of

with multipliers $\lambda_{i}$, $i=1,\dots,k$, and finding the values of $\lambda_{i}$ for which equalities in (3) hold. It is easy to see that (5) is a particular case of (4) with $\lambda_{ij}=\lambda_{i}$ for all $j=1,2,\dots,k$, $j\not=i$, so in what follows we focus on the minimization of (4).

It is not difficult to see that in the particular case when $\theta_{i}=\vartheta_{i}$, $i=1,2\dots,k=K$ the Lagrangian function (4) can be considered Bayesian risk (see, for example, Baum and Veeravalli 1994, among many others) corresponding to the a priori distribution $(\gamma_{1},\dots,\gamma_{k})$ on the set of parameter points $\{\theta_{1},\dots,\theta_{k}\}$, where $\lambda_{ij}/\gamma_{i}$ can be interpreted as conditional loss from accepting $H_{j}$ when $H_{i}$ is true. Thus, the minimization of (4) readily solves the problem of optimal Bayesian tests for $k$ hypotheses.

The well-known modified Kiefer-Weiss problem (see, for example, Lorden 1980) also easily embeds into this scheme by taking $\gamma_{1}=1$, $K=1$, and $\vartheta_{1}$ between the hypothesized values $\theta_{1}$ and $\theta_{2}$, being $k=2$. And this gives rise to a multi-hypothesis version of the Kiefer-Weiss problem, starting from a modified version of it, with $\vartheta_{1},\dots,\vartheta_{k-1}$ such that $\theta_{1}<\vartheta_{1}<\theta_{2}<\vartheta_{2}<\dots<\vartheta_{k-1}<\theta_{k}$ and with some weights $\gamma_{1},\gamma_{2},\dots,\gamma_{k-1}$, adding up to 1, as additional parameters. To our knowledge, there are no known non-asymptotic solutions of the multi-hypothesis Kiefer-Weiss problem, and this could be a basis for one.

Now, let us characterize the tests which minimize the Lagrangian function (4), for a given set of multipliers. It is worth noting that $L(\psi,\phi)$ implicitly depends on the Lagrange multipliers, therefore all the constructions below will also (implicitly) depend on $\lambda_{ij}$, as well as on other elements of problem setting, like $\theta_{i}$ and $\vartheta_{i}$, etc.

First of all, in a very standard way it can be shown that there is a universal decision rule $\phi$ that minimizes $L(\psi,\phi)$ whatever fixed $\psi$ (see Novikov 2009b).

Let us assume that $P_{\theta}$ is absolutely continuous with respect to a $\sigma$-finite measure $\mu$ and denote $f_{\theta}$ its Radon-Nikodym derivative. Also denote $f_{\theta}^{n}=f_{\theta}^{n}(x_{1},\dots,x_{n})=\prod_{i=1}^{n}f_{\theta}(x_{i})$, and let $f^{n}_{{\gamma\vartheta}}=\sum_{i=1}^{K}\gamma_{i}f_{\vartheta_{i}}^{n}$. Define

Let a decision rule $\phi$ be such that

(in the case of equality in (7) $\phi_{n}^{j}$ can be arbitrarily randomized between those $j$ sharing this equality, with the only requirement that $\sum_{j=1}^{k}\phi_{n}^{j}\equiv 1$). It follows from Theorem 3 in Novikov (2009b) that

and we have an optimal stopping problem of minimizing (8) over stopping rules $\psi$.

The problem is first solved in the class of truncated tests, i.e. those not taking more than a finite number $N$ of observations. Let $\mathcal{S}^{N}$ be the set of all such stopping rules that $(1-\psi_{1})\dots(1-\psi_{N})\equiv 0$.

Let us define operator $\mathcal{I}_{n}$ in the following way. For any measurable non-negative $v=v(x_{1},\dots,x_{n})$ let

Now, starting from

define recursively over $n=N,N-1,\dots,2$

Then for any $\psi\in\mathcal{S}^{N}$

and there is an equality in (9) if for all $n=1,2,\dots,N-1$

where $I_{A}$ denotes the indicator function of the event $A$. In this way, stopping rule $\psi$ in (10) minimizes $L(\psi)$ in the class of truncated stopping rules $\mathcal{S}^{N}$. Any $\psi_{n}$ may be arbitrarily randomized between samples $(x_{1},\dots,x_{n})$ for which there is an equality in the inequality under the indicator function in (10). This gives the same value of $L(\psi)$. The details can be found in Novikov (2009b).

The optimal non-truncated tests can be found passing to the limit, as $N\to\infty$, provided that

(see Remark 7 in Novikov 2009b). In the case of i.i.d. observations we are considering in this paper, (11) holds without any additional conditions. The formal proof of this fact can be found in the Appendix.

The construction of the optimal non-truncated test is as follows. First of all, it is easy to see that $V_{n}^{N+1}\leq V_{n}^{N}$, so there exists $V_{n}=\lim_{N\to\infty}V_{n}^{N}$, $n=1,2,\dots$ Then it follows from (9) that

and the right-hand side in (12) is attained if

for all $n=1,2,\dots$ In this way, we obtain tests $\langle\psi,\phi\rangle$ with $\psi$ satisfying (13) and $\phi$ satisfying (7) which minimize the Lagrangian function $L(\psi,\phi)$.

We propose using numerical methods for construction of the truncated tests minimizing the Lagrangian function. For the Bernoulli model, we develop numerical algorithms for this and implement them in the form of a computer program in the R programming language. Having the means for minimizing the Lagrangian function, to obtain optimal sequential tests in the conditional setting (i.e. those minimizing $C_{\gamma,\vartheta}$ under conditions (2)) we need to find Lagrangian multipliers $\lambda_{ij}$, $1\leq i\not=j\leq k$, providing a test (7)-(10) for which equalities in (2) hold. Respectively, the minimization of $C_{\gamma,\vartheta}$ under conditions (3) reduces to finding $\lambda_{i}$, $i=1,\dots,k$ such that for the test in (7)-(10), with $\lambda_{ij}=\lambda_{i}$ for $1\leq j\not=i\leq k$, for which there are all equalities in (3).

In no way can one be sure that such $\lambda_{ij}$ exist for every combination of $\alpha_{ij}$ (not even in the classical case of two hypotheses). On the other hand, every combination of $\lambda_{ij}$ employed in (7)-(10), produces an optimal test $\langle\psi,\phi\rangle$ in the conditional setting, if one takes its error probabilities as $\alpha_{ij}$ in (2) (i.e. $\alpha_{ij}=\alpha_{ij}(\psi,\phi)$) (or, respectively, as $\alpha_{i}$ in (3), that is $\alpha_{i}=\alpha_{i}(\psi,\phi)$).

Having at hand a computer program for the Lagrange minimization, finding the multipliers providing a tolerable level of the error probabilities is a question of some trial-and-error look-ups, because larger values of $\lambda_{ij}$ make $\alpha_{ij}$ smaller, grosso modo. As an alternative, general-purpose computer algorithms of numerical optimization can be used to get as close as possible to the desired values of $\alpha_{ij}$ by moving the input values of $\lambda_{ij}$, for example, the method of Nelder and Mead (1965).

For the non-truncated tests, we propose using approximations by truncated tests. We illustrate all this technique on the particular case of Bernoulli distribution in the subsequent sections.

In this section, we apply the general results of Section 2 to the model of Bernoulli observations. In this way we obtain a complete set of computer algorithms for computing the tests that minimize the Lagrangian function, and their numerical characteristics, in the Bernoulli model. For the determination of the values of the Lagrange multipliers general-purpose computer algorithms will be used.