Bayesian Sequential Joint Detection and Estimation under Multiple Hypotheses
Dominik Reinhard Michael Fauß Abdelhak M. Zoubir
eess.SP cs.IT math.IT stat.AP







Bayesian Sequential Joint Detection and Estimation under Multiple Hypotheses
    Alexey BufetovInstitute for Applied Mathematics, Bonn University, Endenicher Allee 60, 53115 Bonn, Germany. E-mail: alexey.bufetov@gmail.com
Patrik L. FerrariInstitute for Applied Mathematics, Bonn University, Endenicher Allee 60, 53115 Bonn, Germany. E-mail: ferrari@uni-bonn.de

    
============================================================================================================================================================================================================================================================================================





[0]Address correspondence to Dominik Reinhard,
	Signal Processing Group, Technische Universität Darmstadt, Merckstraße 25, 
	64283 Darmstadt, Germany; E-mail: reinhard@spg.tu-darmstadt.de


[1]inputs: #1

[1]initialize: #1


Abstract:

 We consider the problem of jointly testing multiple hypotheses and estimating a random parameter of the underlying distribution. This problem is investigated in a sequential setup under mild assumptions on the underlying random   process.
 The optimal method minimizes the expected number of  samples while ensuring that the average detection/estimation  errors do not exceed a certain level.
 After converting the constrained problem to an unconstrained one, we characterize the general solution by a non-linear Bellman equation, which is parametrized by a set of cost coefficients. A strong connection between the derivatives of the cost function with respect to the coefficients and the detection/estimation errors of the sequential procedure is derived. 
 Based on this fundamental property, we further show that for suitably chosen cost coefficients the solutions of the constrained and the unconstrained problem coincide.
 We present two approaches to finding the optimal coefficients. For the first approach, the final optimization problem is converted into a linear program, whereas the second approach solves it with a projected gradient ascent.
 To illustrate the theoretical results, we consider two problems for which the optimal schemes are designed numerically. Using Monte Carlo simulations, it is validated that the numerical results agree with the theory.
 



Keywords: 
Joint detection and estimation, Multiple hypotheses; Sequential analysis; Stopping time



Subject Classifications: 
	62L10; 62L12; 62L15; 90C05; 93E10.



§ INTRODUCTION

In many applications, hypothesis testing and parameter estimation occur in a coupled way and both are of equal interest. More precisely, one has to decide between two or more hypotheses and, depending on the decision, estimate one or more, possibly random, parameters of the underlying distribution. This problem goes back to Middleton et al., who initially addressed this problem in the late 1960s. In <cit.>, a Bayesian framework was used to find a jointly optimal solution. After extending the problem to multiple hypotheses <cit.>, the interest in the topic declined in the literature, but the topic regained more attention <cit.> since 2000. 

Joint detection and estimation is of particular interest in signal processing and communications. For example, in cognitive radio, the secondary user has to detect the primary user and estimate the possible interference such that dynamic spectrum access can be performed <cit.>. In a general communication setup, one is interested in detecting the presence of a signal and estimating the channel<cit.>.
There exist many more applications, such as radar<cit.>, speech processing <cit.>, change point detection and estimation of the time of change <cit.>, optical communications <cit.>, detection and estimation of objects from images <cit.> and biomedical engineering <cit.>.

In applications such as radar, where the aim is to detect a single target and to estimate, for example, the velocity of the target, formulating the joint detection and estimation problem with two hypotheses is sufficient.
Nevertheless, there exist a variety of applications where we have to decide between multiple hypotheses and to estimate simultaneously one or more parameters.
In communications, for example, one wants to decode a received signal and estimate some parameters of the signal model <cit.>, e.g., the unit sample response of the channel or the noise power.
If the alphabet consists of more than two symbols, testing multiple hypotheses is required.
Another example is state estimation for smart grids, where topological uncertainties can be modeled using multiple hypotheses <cit.>. 

Constraints on the latency of a system arise in many signal processing applications.
Therefore, it is required to perform inference as fast as possible and ensure at the same time its quality.
This leads to the field of sequential analysis that was pioneered by Wald in the late 1940 <cit.>.
Sequential methods stop sensing new data once one is confident enough about the phenomenon of interest.
Compared to their fixed-sample-size counterparts, these methods require on average significantly less samples while they reach the same inference quality.
An overview on sequential detection and estimation can be found in <cit.> and <cit.>, respectively.




Combining the key ideas of sequential analysis with those of joint detection and estimation leads to a powerful framework which is applicable for a wide range of signal processing tasks. In this framework, the average number of used samples should be minimized, while the detection and estimation errors are controlled.




 §.§ State of the Art

The design of optimal procedures for the problem of sequential joint detection and estimation is an important yet little studied topic.
There exist well established sequential tests in the literature that include an estimation step, such as the generalized SPRT or the adaptive SPRT.
For those tests, only the outcome of the hypothesis test is of interest and hence, it is not possible to control the quality of the estimate.
Details about the generalized SPRT and the adaptive SPRT can be found in <cit.>.

In <cit.>, the authors investigated the sequential update of sufficient statistics for simultaneous sequential detection and estimation. However, the authors did not provide an optimal solution for the joint detection and estimation problem.  
The problem of sequential detection and state estimation for Markov processes was investigated in <cit.> .
Although the quality of the estimator, in that case the probability of a correct estimation, is of primary interest, it does not affect the number of used samples.
Later on, a test that incorporates both uncertainties, the one about the true hypothesis and the true states, into the number of used samples was designed for a similar problem <cit.>.

The design of optimal procedures for the problem of sequential joint detection and estimation was initially treated by Yılmaz et al. In <cit.>, the aim is to decide between two hypotheses and, if the null is rejected, to estimate a random parameter. The approach in <cit.> was later extended to multiple hypotheses <cit.> and applied to joint spectrum sensing and channel estimation <cit.>. The objective in the work by Yılmaz et al. <cit.> is to design a scheme that minimizes the number of used samples for every set of observations while keeping a combined detection/estimation cost function below a certain level. The cost function is a linear combination of detection and estimation errors with combination weights set by the designer. Although the detection and estimation errors can be balanced by varying the combination weights, the main drawback of that approach is the difficulty in designing  procedures with predefined performance in terms of error probabilities and estimation errors.

In <cit.>, we investigated the problem of sequential joint detection and estimation under distributional uncertainties. In <cit.>, the problem of sequential joint signal detection and SNR estimation was addressed. In that work, the SNR was modeled as a deterministic unknown parameter and the detection and estimation errors were kept below predefined limits.
More recently, we developed a Bayesian framework for two hypotheses.
In that work, the average number of used samples is minimized while the detection and estimation errors are kept below predefined levels<cit.>. Contrary to the work by Yılmaz et al. <cit.> which only allows a balance of detection and estimation errors, the framework in <cit.> comes with an approach to choose the cost coefficients such that constraints on the detection and estimation errors are fulfilled and such that the resulting scheme uses on average as few samples as possible.
The framework was later applied to joint signal detection and SNR estimation <cit.>.
Sequential joint detection and estimation in distributed sensor networks was addressed in <cit.> and <cit.>.








 §.§ Contributions and Outline

In this work, we show how the framework in <cit.> can be recast to the general case of multiple hypotheses and as a result we unify the theory so that the work in <cit.> becomes a special case of this work.
The aim is to design a strictly optimal method, i.e., a method that fulfills constraints on the detection and estimation errors while it uses on average a minimum number of samples.
To the best of our knowledge, the design and analysis of strictly optimal procedures for sequential joint detection and estimation for the multiple hypotheses case has not been addressed in the open literature so far.
Although this problem might look similar to the one in <cit.>, the multiple hypotheses case differs significantly from its binary counterpart.
These differences are captured in the cost for stopping the scheme and, as a consequence, also in the optimal cost function that characterizes the procedure.
To design the procedure, the coefficients that parameterize the optimal cost function have to be chosen such that the constraints on the error probabilities and the MSE are fulfilled.
To this end, a connection between the performance measures, i.e., the detection and estimation errors, and the partial derivatives of the cost function with respect to the corresponding coefficients is exploited as in <cit.>.
However, due to the structural differences of the optimal cost function in the binary and the multiple hypotheses scenario, a connection between the error probabilities and the partial derivatives of the cost function is a non-obvious one.
The proof of this property is given in this work. Although similar tools as in <cit.> are used, the two proofs differ technically. 
Based on this connection, the original design problem can be, as its binary counterpart, converted to a linear program.
In addition, we provide a second design algorithm that is based on a projected gradient ascent which was not a subject in <cit.>.

Note that a small amount of some results in this paper were presented at the 2020 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)<cit.>.

This work is structured as follows: in <ref>, we summarize the underlying assumptions and provide a formal problem formulation. <ref> describes the conversion of this problem to an optimal stopping problem. The solution of this problem, which can be obtained by dynamic programming, is characterized by a set of non-linear Bellman equations. The properties of the resulting cost functions are discussed in <ref>. In <ref>, we present two approaches for selecting the  coefficients, which parametrize the cost functions, such that predefined detection and estimation performances are achieved and the resulting scheme is of minimum average run-length. To illustrate the proposed approaches, <ref> provides two numerical examples. 
 

§ PROBLEM FORMULATION

Let = (X_n)_n ≥ 1 be a sequence of random variables, which is generated under one out of M different hypotheses _m, m=1,…,M. Under each hypothesis, the distribution of the sequence  is parametrized by a single scalar and real random parameter _m with known distribution.
Although the problem can be formulated with multiple parameters under each hypothesis, for the sake of simplicity we focus on the single parameter case. The work on the multiple parameter case will appear elsewhere.
The tuple (, , ) is defined on the probability space (××, ℬ() ⊗ℬ(ℝ) ⊗ℱ_ , P) in the metric state space (××, ℰ_X×ℰ_×ℰ_). The density of P with respect to a
suitable reference measure, for example, the standard Lebesgue measure when  is Euclidean, is denoted by p.
Hence, the M composite hypotheses can be written as

    _m: _m ∼ p(𝐱_m, _m) , _m∼ p(_m_m) ,

where m=1,…,M and 𝐱 are the observations of the random sequence 𝒳.
Note that the underlying hypothesis can, in the most general case, affect both, the distribution of the data given the parameter, as well as the distribution of the parameter.
Since the M hypotheses are composite, we do not only want to decide for the true hypothesis, but we are also interested in estimating the underlying parameter _m.
Hence, the problem is one of joint detection and estimation. Moreover, our aim is to solve this problem in a sequential setup, i.e., we observe the sequence  sample by sample until we are confident enough about the true hypothesis and the underlying parameter. Thus, the problem becomes one of sequential joint detection and estimation. 
The decision maker is restricted to take at most N samples.
That is, irrespective of the certainty about the true hypothesis and the true parameter, the sequential scheme has to stop sampling and to infer the quantities of interest.
For a brief discussion on how to choose N in practice, see <ref>.


Before going into the details of the optimization problem, we first summarize the assumptions for the proposed framework. Then, some fundamentals of sequential joint detection and estimation are provided, followed by a discussion of the performance metrics.



 §.§ Assumptions and Notation



  * The true hypothesis and the true parameter do not change during the observation period.

  * The M hypotheses are mutually exclusive.
  *  There exists a sufficient statistic n() in a metric state space (, ℰ_) such that for all n=1,…,N, m=1,…,M,

    p(n+1, _m, _m ) = p(n+1, _m, _m | n)

with some initial statistic 0 and a transition kernel

    n+1(n+1) = ξ(n,n+1) =: (n+1) .


  * The second order moments of the random parameters _m conditioned on the hypotheses _m exist and are finite:

    [_m^2_m] < ∞ ,  ∀ m=1,…,M



For the sake of compactness, the integration domain is not indicated explicitly in integrals that are taken over the whole domain, e.g.,

    [X] = ∫_ xp(x) x ≡∫ xp(x) x .

Moreover, the dependency of functions, such as estimators, on the data or the sufficient statistic is dropped for simplicity. The dependency should be clear from the context.
The indicator function of the event 𝒜 is denoted by 𝒜.



 §.§ Fundamentals of Sequential Joint Detection and Estimation

Since the sequence  is observed sample by sample, the state space of the observations grows with every time instant. To keep the problem tractable, the observations n are replaced by the sufficient statistic n as stated in <ref>, which serves as a low-dimensional representation of the data. Moreover, its state space is independent of the time instant n.


For the joint detection and estimation problem, a decision rule, as well as a set of estimators (one under each hypothesis) have to be found. The decision rule maps the sufficient statistic to a decision in favor of one particular hypothesis.
The estimators, on the other hand, map the sufficient statistic onto a point in the state space of the parameter under each hypothesis. Mathematically, the decision rule and the estimators are defined as

    _n    : ↦Δ_M = {1,…,M} ,
    m,n   : ↦ ,    m=1,…,M .

The decision rule and the estimators depend on the number of samples n, which is indicated by the subscript.
In a sequential framework, the number of used samples is not given a priori. Thus, a stopping rule needs to be found for deciding whether to stop sampling or not. Mathematically, the stopping rule is defined as

    _n    : ↦Δ = {0,1} ,

where 0 and 1 stand for continue sampling and stop sampling, respectively. Hence, the run-length of the sequential scheme can be defined as

    τ = min {n≥ 1: _n = 1} .

Note that, since the sopping rule evaluates the sufficient statistic, which is a random variable, the run-length is also a random variable.
Contrary to, e.g., <cit.>, where the stopping and decision rules are probabilities instead of hard decisions, there is no need for randomization in this work. The reason will become clear later.
The collection of the decision rule, the estimators and the stopping rule is referred to subsequently as policy, and is defined as

    =  .

The set of all feasible policies is given by

    = {: ↦Δ^N ×Δ_M^N ×^MN} .



 §.§ Performance Measures

One common performance measure for sequential inference is the average run-length, i.e., the number of samples which are used on average.
Besides that, there exist different performance measures when dealing with multiple hypotheses, e.g., the probability of accepting _i when _m is true P({_τ = i}_m), i≠ m, or the probability of falsely rejecting hypothesis _m P({_τ≠ m }_m).
In this work, we focus on the latter although an extension to the former should not be too difficult.
As for the estimation performance, there exist many measures such as the MSE or the MAE. In this work, the MSE is used to quantify the estimation performance.
More formally, the error metrics can be written as

    nm(n)     =  [_τ≠ m_m, n, τ≥ n] , 
    nm(n)     =  [_τ= m(_m - m,τ)^2_m, n, τ≥ n] ,

where

    {_τ=m} := {n∈: _τ(n) = m} .

This short hand notation for subspaces of the state space is used throughout the paper.

The quantity nm(n) denotes the probability that hypothesis _m is rejected given that _m is true and that the scheme is in state n at time n.
Similarly, nm(n) is the MSE under the condition that the scheme is in state n at time n. The estimation error is set to zero if the scheme decides erroneously for a hypothesis.

Alternatively, we can express the performance measures in <ref> and <ref> recursively. These definitions, which follow from the Chapman-Kolmogorov backward equations <cit.>, are given by

    nm(n)     =  _n_n≠ m  + (1-_n) [n+1m(n+1)_m, n] ,
    Nm(N)     =  _N_N≠ m , 
    nm(n)     =  _n_n= m[(_m - m,n)^2_m,n]  + (1-_n)[n+1m(n+1)_m, n] , 
    Nm(N)     =  _N_N= m[(_m - m,N)^2_m,N] .

The importance of these recursive definitions will become clearer later. For the sake of compactness, the following short-hand notations will be used throughout the paper

    m := 0m(0) ,  m := 0m(0) .





 §.§ Formulation as an Optimization Problem

The aim is to design a sequential method which is of minimum expected run-length while the detection and the estimation errors do not exceed certain levels.
As mentioned earlier, the decision maker is allowed to observe at most N samples, i.e., the sequential scheme is truncated. The maximum number of samples might be determined by the application, e.g., due to latency constraints.
If the choice is left to the designer, it should be chosen as large as possible given the computational resources.
However, if N is chosen too small, i.e., the constraints on the detection and/or the estimation errors cannot be fulfilled, the design problem does not have a feasible solution and the algorithm(s) presented in <ref> will not converge.

Formally, the problem can be written as the following optimization

    min_∈ [τ] , _N = 1 ,
    subject to     m≤ ,  m≤ ,    m=1,…,M ,

where the upper limits on the detection and estimation errors are denoted by  and , respectively.
Instead of solving <ref> directly, we first consider the following unconstrained problem


    min_∈ {[τ] + ∑_m=1^M( _m p(_m)m + _mp(_m)m)} ,

where _m and _m, m=1,…,M, are some non-negative and finite cost coefficients. Although the prior probabilities p(_m) are not subject to optimization, they are included to simplify notation. It is shown in <ref> that for suitable chosen _m and _m, the solutions of <ref> and <ref> coincide.
 

§ REDUCTION TO AN OPTIMAL STOPPING PROBLEM

In order to solve the problem in <ref>, it first has to be converted to an optimal stopping problem. Optimal stopping theory is a framework which is widely used when designing sequential schemes. In that framework, a stopping rule has to be found which trades-off the expected run-length and the inference quality.  More details about optimal stopping and its applications are given in <cit.>.
To derive an optimal stopping problem, we proceed as in <cit.> and optimize <ref> first with respect to the decision rule and then with respect to the estimators. To do so, the two summands in <ref>, the expected run-length and the combined detection and estimation errors have to be reformulated. 
The expected run-length can be expressed as

    [τ]     = ∑_n=1^N ∑_m=1^M ∫∫ n  p(_m, _m, )_m
        = ∑_n=1^N ∫ n  p() ,

where we used the short-hand notation

    := _n∏_k=0^n-1(1-_k) .


Similarly, the weighted sum of detection and estimation errors in <ref> can be simplified to

    ∑_m=1^M( _m p(_m)m + _mp(_m)m) =    ∑_n=0^N ∫∑_m=1^M  _n=m D_m,n(n()) p() .

The cost for stopping at time n and deciding in favor of _m, D_m,n(n), is given by

    D_m,n(n)     = _m p(_mn)∫(_m -  m,n)^2p(_m_m,n)_m + ∑_i=1,i≠ m^M _i p(_in) .

Hence, the overall objective becomes

    min_∈ ∑_n=0^N ∫(n+∑_m=1^M _n=m D_m,n)p() .
In order to solve <ref>, it is first minimized with respect to the decision rule and then with respect to the estimators.
Since D_m,n(n) is non-negative for all n=0,…,N, m=1,…,M, and all n∈, it holds that <cit.>

    ∫∑_m=1^M   _n=m D_m,np()≥∫min_mD_m,n p() ,

where equality holds if and only if

    _n(n) ∈{m: D_m,n(n) = min_1≤ i ≤ M D_i,n(n)} .
In general, more than one hypothesis may fulfill <ref> which then calls for randomizing _n. In this work, there is no need for randomization. The reason is the same as that stated in <cit.>.
After we derived the optimal decision rule, the optimal estimators have to be found. Since the optimal sequential estimator is independent of the stopping time <cit.> and each D_m,n only depends on one estimator, we can minimize all D_m,n with respect to each estimator m,n separately.
The estimator which minimizes the MSE, i.e., the MMSE estimator, is 
    m,n^⋆ = [_m_m,n]
and the MMSE is given by the posterior variance <cit.>. Hence, the cost for stopping and deciding in favor of _m at time n becomes

    D_m,n^⋆(n)     = _m p(_mn) [_m_m, n]  + ∑_i=1,i≠ m^M _i p(_in) .
The optimal stopping rule can now be found by solving the following optimization problem

    min_{_n}_n=0^N∑_n=0^N [(n + g(n))] ,

where

    g(n) = min{D^⋆_1,n(n) , … , D^⋆_M,n(n)}

is the cost for stopping at time n. The solution of the optimization problem is fixed in the following theorem.
 The solution of <ref> is characterized by the non-linear Bellman equations
 
    ρ_n(n)     = min{g(n), d_n(n)}   n<N  ,
    ρ_N(N)     = g(N) ,

 with g(n) defined in <ref> and the cost for continuing is
 
    d_n(n) = 1 + ∫ρ_n+1(())p(n) .


Since the proof of <ref> does not differ from the ones in the literature, we refer to, e.g., <cit.>,<cit.>.
With a change in measure, <ref> becomes

    d_n(n) = 1 + ∫ρ_n+1 ,
where

    (B) := P({∈: () ∈ B}n)
for all elements B of the Borel σ-algebra on . The probability measure ^m, which will be needed later, is defined analogously, but with P(·n) replaced by P(·_m,n).

 Let _m and _m be finite for all m=1,…,M, then ρ_n+1 is -integrable for all n in  and all 0≤ n<N.

A proof of <ref> is given in <ref>.
The optimal policy is summarized in the following corollary.

 The optimal policy which solves <ref> is 
 
    =  ,

with _n defined in <ref> and m,n^⋆ defined in <ref>. The optimal stopping rule _n is given by

    _n(n) = g(n) = ρ_n(n) .


For the optimal policy stated in <ref>, the stopping region of the scheme, its complement and its boundary are

    = {n∈: g(n) < d_n(n)}
        = {n∈: g(n) = d_n(n)}
        = {n∈: g(n) > d_n(n)}

for n<N.
Moreover, let

    nm = ∪{n∈: _n(n) = m}
    nm̅ = ∪{n∈: _n(n) ≠ m}

denote the regions in which the procedure stops and accepts/rejects hypothesis _m.
 

§ PROPERTIES OF THE COST FUNCTION

In this section, we present, similarly to <cit.>, the fundamental properties of the cost functions obtained in the last section. However, as mentioned before, these properties do not follow directly from their two-hypotheses counterpart.
These properties are used later to get the optimal cost coefficients such that the solution of <ref> also solves <ref>. In order to simplify the upcoming derivations, it is important to show that the boundary of the stopping region is a P-null set. This is stated in the following lemma.

 If the posterior probabilities p(_mn), m=1,…,M, are continuous random variables with respect to n, it holds that
 
    () = 0 ,  ∀n∈,  0≤ n < N ,

i.e., the boundary of the stopping region  is a P-null set.

The proof of <ref> is laid down in the supplement <ref>. <ref> implies that the cost minimizing stopping rule is unambiguous so that there is no need for randomization. Before the main properties of the cost functions are presented, the following short hand notations are introduced

    z_n^m     := p(n_m)/p(n) ,
    {∈[n+1]}    := {∈: ()∈[n+1]} .


Let ρ_n,_m^' and ρ_n,_m^' denote the derivatives of ρ_n with respect to _m and _m for m=1,…,M, respectively. For n<N, it holds that
 
    ρ_n,_m^'(n)     =  nm̅ p(_mn)  + ( p(_m)z_n^m^m(n+1m̅) + ∫_{∈[n+1]}ρ^'_n+1,_m)

 and
 
    ρ_n,_m^'(n)     = nm p(_m) z_n^m [_m_m, n] + r_n^m

 with r_n^m being recursively defined via
 
    r^m_n     = p(_m) z_n^m ∫_{∈n+1m}[_m_m,()]p(n,_m) + ∫_{∈[n+1]}ρ^'_n+1,_m .

For n=N, it further holds that

    ρ_N,_m^'(N)     = Nm̅p(_m)z_N^m 
    ρ_N,_m^'(N)     = Nm p(_m) z_N^m [_m_m, N] .


The connection of the derivatives with respect to _m do not follow from <cit.> straightforwardly. Therefore, the proof is given in <ref>. For the derivative with respect to _m, we refer to <cit.>.
<ref> is an intermediate result, which is required to prove the more important result stated in <ref>.
Next, we provide a connection between the derivatives of the cost functions stated in the previous theorem and the performance measures stated in <ref>. This connection forms the basis to obtain the optimal cost coefficients, as presented in <ref>.

Let ρ_n,_m^' and ρ_n,_m^' be as defined in <ref>. Then, using the optimal policy given in <ref>, it holds that 
 
    ρ_n,_m^'(n)     = p(_m)z_n^m nm(n) , m=1,…,M ,
    ρ_n,_m^'(n)     = p(_m)z_n^m nm(n) , m=1,…,M ,

 and in particular
 
    ρ_0,_m^'(0)     = p(_m)0m(0) , m=1,…,M ,
    ρ_0,_m^'(0)     = p(_m)0m(0) , m=1,…,M .


As the connection for the derivatives with respect to _m and the detection errors do not follow from its two-hypotheses counterpart, the proof is outlined in  <ref>. For the derivatives with respect to _m, we refer to <cit.>.
 

§ CHOICE OF THE OPTIMAL COST COEFFICIENTS

This section extends Section 5 in <cit.> to multiple hypotheses.
Proofs are given for results that do not follow directly from <cit.>.

Often, the question how to choose the coefficients of a cost function stays untouched and the choice is left to the designer <cit.>. However, in this work, our aim is to design a sequential scheme which is of minimum average run-length and which fulfills predefined constraints on the detection and estimation errors. Since all three performance measures, average run-length, error probabilities and MSE, are of different numerical range, it is rather impossible to choose the cost coefficients by hand.
In order to automatically select the correct cost coefficients such that a predefined performance is achieved, the results stated in <ref> are exploited.

For the sake of a compact notation, let λ=(λ_1,…,λ_M) and μ=(μ_1,…,μ_M). To obtain the set of optimal cost coefficients, we consider the following maximization problem, which is in fact the Lagrangian dual of <ref>,

    max_≥0,≥0     L_[],[](,) ,

where ≥0, ≥0 have to be read element-wise and 

    L_[],[](,) = ρ_0(0) - ∑_m=1^M p(_m)(_m + _m) .

As it is not trivial to see that strong duality holds and therefore the solutions of <ref> and <ref> coincide, it is fixed in the following theorem.

Let  be the solution of <ref>, let _[],[] and _[],[]  be solutions of <ref> and let  be the policy parametrized by _[],[], _[],[]. Then, it holds that

    =  , 
    
     L_[],[](_[],[], _[],[])     = [τ=] .

That is, a solution of <ref> also solves <ref>. Moreover, the optimal objective of <ref> is the expected run-length.

A proof of <ref> is outlined in <ref>.

Hence, by using <ref> and <ref> the original problem given in <ref> is equivalent to

    max_≥0,≥0    {ρ_0(0) -∑_m=1^M p(_m)(_m + _m)} , 
    s.t.    ρ_n(n) = min{g(n), d_n(n)} ,   n<N , 
       ρ_N(N) = g(N) .

In what follows, we present two approaches to solve <ref>. The first one uses linear programming to obtain the optimal cost coefficients and the cost functions, whereas the second one uses a projected gradient ascent.



 §.§ Linear Programming

The first approach solves <ref> by linear programming. To this end, we proceed similarly to <cit.> and relax the equality constraints in <ref> to multiple inequality constraints and add the cost functions to the set of free variables, i.e.,

    max_≥0,≥0
    ρ_n∈ℒ    {ρ_0(0) - ∑_m=1^Mp(_m)(_m + _m)} ,
    s.t.     ρ_n ≤ D^⋆_m,n ,  m=1,…,M,  n=0,…,N , 
       ρ_n ≤ 1+∫ρ_n+1 ,   n=0,…,N-1 ,

where ℒ is the set of all -integrable functions on .

 Problem <ref> is equivalent to Problem <ref>.

For a proof see <cit.>.


 §.§ Projected Gradient Ascent

Besides the previously presented LP, we present a second approach that is not as sensitive to the size of the design problem as the LP. This approach uses a projected gradient ascent to obtain the optimal cost coefficients and is summarized in <ref>.

At each iteration, the cost functions, as stated in <ref>, are calculated for a given set of cost coefficients. Next, given the policy induced by these cost functions, the gradient of the objective in <ref> has to be obtained, which is given by

    ∇_ L_[],[](, )     = [p(_1)(1 -[1]),…, p(_M)(M-[M]) ] ,
    ∇_ L_[],[](, )     = [p(_1)(1 -[1]),…, p(_M)(M-[M]) ] .

This gradient can, e.g., be calculated based on the definitions of the performance measures as stated in eq:perfMeasuresRecFirsteq:perfMeasuresRecLast.
To update the cost coefficients, the old coefficients are shifted in the direction of the gradient and then projected onto the set of feasible coefficients, i.e., the set of non-negative reals.
To control the convergence speed of the algorithm, the gradient is scaled by a factor .
These steps are repeated until the solution converges to an optimum. Recall from <ref> that the optimality criteria for the cost coefficients are

    _m(m - )     = 0 ,  m=1,…,M ,
    _m(m - )     = 0 ,  m=1,…,M .

Hence, the procedure has to be repeated until <ref> holds (approximately).

Since the calculation of the performance measures based on their recursive definition can become numerically unstable and hence lead to an inaccurate gradient, a modification of the aforementioned method can be used. Similarly to <cit.>, the detection and estimation errors in <ref> and <ref> can be replaced by their Monte Carlo estimates.


As Problem <ref> is convex, the projected gradient ascent converges to a global optimum irrespective of the choice of the starting point. Nevertheless, the choice of the starting point is crucial for the convergence speed. To have a fast convergence, we suggest to solve <ref> on a coarse grid if possible and then run <ref> on a finer grid to obtain the optimal cost coefficients.






 §.§ Discussion

Although the algorithms presented in this section are in general problem independent, their tractability highly depends on the size of the problem.
With increasing dimensionality of the sufficient statistic, the discretization of its space and hence, the calculation of the optimal policy becomes infeasible.
However, even if the observations are high-dimensional, it might still be possible to find a low-dimensional sufficient statistic.
An example for such an application could be a sensor network in which all sensors send their observations to a fusion center that performs the inference.

Besides the dimensionality of the sufficient statistic, also the dimensionality of the observations may play a role in the implementation of the algorithms.
The cost for continuing, i.e., d_n(n) = 1 + [ρ_n+1(n+1)n], can be 
calculated numerically in a straightforward manner for the case of scalar observations. For more dimensional observations, approximating the integral by Markov Chain Monte Carlo integration, as, e.g., in <cit.>, might be more suitable.

Nevertheless, to overcome this shortcoming, more sophisticated design approaches that do not require a discretization of the state space have to be found.
We have ongoing work on machine learning techniques to learn the policy from observed trajectories rather than calculating the cost function on a discretized state space.

 

§ NUMERICAL RESULTS

In this section, we provide two numerical examples to illustrate and validate the proposed approach. First, a simple example is presented to illustrate the basic properties of the optimal method. The second example is more complex and shows how to apply the proposed method to real-life applications.

In order to solve the linear program in <ref>, the continuous spaces are first discretized. The discrete linear program is then solved by the Gurobi optimizer <cit.> which is called via the MATLAB cvx interface <cit.>.
We add a regularization term to the objective in <ref> to ensure numerical stability.
See, e.g., <cit.> for details. To validate the performance of the designed schemes, a Monte Carlo simulation is performed for both examples. 



 §.§ Benchmarking Method

In order to compare the proposed approach with existing methods, we choose a two-step procedure as benchmarking method, namely, a standard sequential detector for multiple hypotheses followed by an MMSE estimator.
Although there exist different sequential detectors for the multiple hypotheses case, we resort to the MSPRT<cit.> due to its easy implementation and favorable theoretical properties. Alternative sequential detectors for multiple hypotheses can be found in, e.g., <cit.>.
This two-step procedure is not optimal for the joint detection and estimation problem, but the MSPRT is asymptotically optimal for the detection part and the MMSE estimator is the optimal estimator with respect to the MSE.
For the MSPRT, the pair-wise log-likelihood ratios for hypotheses _m and _j are used and are defined as

    _mj(n) = log(p(n_m)/p(n_j)) ,    
    m,j =1,…,M ,   m  ≠ j .

In general, the stopping and decision rules of the MSPRT are given by<cit.>:

    _n^MSPRT    =
    					  1    ∃ m: _mj≥ A_mj, ∀ j∈{1,…,M}∖ m ,
    
                                              0    else ,
    
    _n^MSPRT    = {
    					  m: _mj≥ A_mj, ∀ j∈{1,…,M}∖ m
                                           } .
The thresholds A_mj, which are used in the stopping and decision rules, have now to be set such that the target error probabilities are met.
To keep the probabilities of falsely rejecting hypothesis _m under a certain level , m=1,…,M, the thresholds have to be calculated as <cit.>

    A_mj = A_m≈log(M/) .

Since the presented optimal sequential scheme is a truncated one, we use a truncated two-step procedure for the sake of fair comparison. Hence, the stopping and decision rules at the truncation point are given by

    _N^MSPRT    = 1   and  _N^MSPRT =  _m ∑_j=1,j≠ m^M _mj .





 §.§ Shift-in-Mean

The first numerical example is used to show the basic properties of the optimal sequential scheme. Here, we consider three different hypotheses with a Gaussian likelihood. Under each hypothesis, the mean follows a different distribution, whereas the variances of the likelihood are equal. The three different hypotheses are given by

    _1:     x_nμ_1 iid∼μ_1σ^2 , -μ_1 + 1.3 ∼(1.7,1) ,
    _2:     x_nμ_2 iid∼μ_2σ^2 , - + 1.3μ_2 ∼(-1,1) ,
    _3:     x_nμ_3 iid∼μ_3σ^2 , -μ_3 - 1.3 ∼(1.7,1) ,

where μσ^2 is the normal distribution with mean μ and variance σ^2, (l,u) is the uniform distribution on the interval [l,u) and (a,b) is the Gamma distribution with shape and scale parameters a and b, respectively.
All three hypotheses have equal prior probabilities and the variance is set to σ^2=4.
The aim is to design an optimal sequential scheme to simultaneously test the three hypotheses and to estimate the mean. The optimal scheme should not use more than 100 samples. The constraints on the detection errors as well as on the estimation accuracy are summarized in the second column of <ref>.
In order to design the optimal scheme, a sufficient statistic in the sense of <ref> has to be found. Let

    x̅_n    = 1/n∑_k=1^n x_k 
    s̅_n^2    = 1/n∑_k=1^n k^2 - x̅_n^2

denote the sample mean and the sample variance, respectively. Then, the likelihood can be written as

    p(   _m,μ_m) = 
        =     (2πσ^2)^-n/2exp(-ns̅_n^2 + n(x̅_n - μ_m)^2/2σ^2)
        = (2πσ^2)^-n/2exp(-ns̅_n^2/2σ^2)exp( - n(x̅_n - μ_m)^2/2σ^2) .

Since the variance is known, the relation between the data and the random mean is completely described by x̅_n, i.e.,

    p(_m,μ_m)    ∝exp( - n(x̅_n - μ_m)^2/2σ^2)
and hence, x̅_n is used as a sufficient statistic in the sense of <ref>. For the likelihood of the sufficient statistic it holds that

    p(n_m,μ_m) ∝exp(-n(x̅_n - μ_m)^2/2σ^2) ∝μ_mσ^2/n .
Note that the likelihood is continuous in the sufficient statistic as well as in the random parameter μ_m, which itself follows a continuous distribution. Therefore, the posterior probabilities p(_mn) are continuous random variables with respect to n and, hence, the boundary of the stopping region is a P-null set according to <ref>.
The discretization of the continuous spaces is summarized in <ref>. 

The posterior probabilities of the hypotheses p(_mn) and the posterior variances [μ_m_m,n] are calculated by numerical integration.
The procedure is then designed by solving the LP in <ref>.

The resulting coefficients are given in <ref> and the optimal policy parametrized by these coefficients as well as the boundaries of the MSPRT  are shown in <ref>. One can see that for small n, the optimal policy is mainly influenced by the detection constraint. For small n, a large value of x̅_n results in a high certainty of _3 but, on the other hand, also in a high estimation error. Hence, the optimal scheme continues sampling in this case and stops only for n>35. Contrary to this, the MSPRT stops directly for large x̅_n, even if n is small, and decides in favor of _3. With increasing n, the optimal policy and the one of the MSPRT become very similar. Nevertheless, the MSPRT has a much broader corridor for continuing than the optimal procedure.


To validate the performance of the optimal and the two-step procedure, a Monte Carlo simulation with 10^6 runs is performed. The results are summarized in <ref> and <ref>.
For the optimal scheme, the constraints are, within the range of Monte Carlo uncertainty, fulfilled with equality and the expected run-length obtained via <ref> is almost equal to the empirical average run-length. For the two-step procedure, the empirical errors under _2 are close to the constraints, but the empirical average run-length of the two-step procedure under _2 is larger than the one of the optimal scheme.
Moreover, under _1 and _3, the average run-length and the detection errors of the two-step procedure are much smaller than the ones of the optimal scheme. 
Though the two-step procedure has smaller average run-lengths and error probabilities under _1 and _3, the empirical MSE is 4 and 8 times as large as the constraint, respectively.

 




 §.§ Joint 4-ASK Decoding and Noise Power Estimation

ASK
This numerical example, which was partly presented in <cit.>, gives an example for how to apply  the proposed framework to a real-world problem.
We consider a 4-ASK (4-ASK) symbol to be transmitted over an additive white Gaussian noise channel with random noise power. At the receiver side, we want to jointly decode the transmitted symbol and estimate the noise power.
Using a linear model, the received signal is given by

    x_n = A + w_n ,

where A∈{A_1,A_2,A_3,A_4} denotes the ASK symbol and w_n, n=1,…,N, is the additive white Gaussian noise process.
According to assumption <ref>, the transmitted symbol A does not change during the observation period. 
The distribution of the noise power σ^2 follows an inverse Gamma distribution itself with known hyperparameters. Here, the signal decoding is a hypothesis test. Hence, the four different hypotheses can be written as

    _m:    σ^2 iid∼A_mσ^2 , σ^2 ∼(a,b) ,

where m=1 …,4, and (a,b) is the inverse Gamma distribution with shape and scale parameters a and b, respectively. The pdf of the inverse Gamma distribution is given by <cit.>

    p(σ^2) = b^a/Γ(a)(σ^2)^-a-1e^-b/σ^2 ,

where Γ(·) denotes the Gamma function.
First, a sufficient statistic in the sense of <ref> has to be found. As shown in <ref>, the likelihood is completely determined by x̅_n as defined in <ref> and s̅^2_n as defined in <ref>.
Hence, the sufficient statistic n=[x̅_n, s̅_n^2] is used in the sequel. 
Since the inverse Gamma distribution is a conjugate prior for the variance of a Gaussian distribution, we can provide analytical expressions for all posterior quantities. First, the posterior distribution of the variance under _m follows itself an inverse Gamma distribution <cit.>, i.e., 
    σ^2_m,n∼(, ) ,
with the posterior parameters

    =  + n/2 ,
        =  + 0.5 ∑_k=1^n (x_k-A_m)^2 
        =  + 0.5 n(s̅_n^2 + (x̅_n - A_m)^2) .

The posterior mean and the posterior variance under _m are then given by <cit.>

    [σ^2_m, n]     = /-1 ,    for >1 ,
    [σ^2_m, n]     = ()^2/(-1)^2(-2) ,  for >2 .

Moreover, an analytical expression for the posterior probabilities of the hypotheses p(_mn) as well as for the posterior predictive p(n) can be expressed as

    p(_mn)     =  K (2π)^-n/2 p(_m) ^/()^Γ()/Γ() , 
    
     p(n)     = K (2π)^-n+1/2∑_m=1^M p(_m) ^/()^Γ()/Γ() ,
where the parameters of the posterior predictive ,  and the normalization constant K are given by

    =  + 0.5 , 
        =  + 0.5(-A_m)^2 , 
    
     K     =(∑_m=1^M p(_mn))^-1 .
A detailed derivation of <ref> and <ref> is laid down in <ref> and <ref>, respectively. Since p(_mn) are continuous random variables with respect to n, the boundary  is a P-null set according to <ref>.

The aim is to design an optimal sequential scheme, which uses at most N=50 samples while the detection and estimation errors are respectively constrained to be below 0.05 and 0.15 under all four hypotheses. The ASK symbols were set to A_m∈{-2,-1,1,2} and the parameters of the noise distribution are given by a=2.1 and b=0.9.
All hypotheses have the same prior probabilities. 

In order to design the optimal scheme, we first use the LP approach on a coarse grid to get an initial set of cost coefficients. These cost coefficients are then used as initial values for the projected gradient ascent which we run on a finer grid.
To reduce the computational load, we exploit the symmetry of the problem at hand. This means that we solve both optimization problems only for _m, _m, m∈{1,2} and ρ_n(n), x̅_n≤0 and complete the missing values once the scheme is designed.


As mentioned before, the LP approach is solved on a coarse grid, whereas a finer grid is used for the gradient ascent. The discretization used for the two algorithms is summarized in <ref>. For the gradient ascent, the gradients are estimated by Monte Carlo simulations with 10^6 runs. The scaling factor for the gradient is set to γ=1000 to speed up convergence.
The stopping criterion for <ref> is

    _m = 0     |m - |   ≤ 10^-3 ,  m=1,…,M ,
    _m = 0     |m - |   ≤ 5·10^-3 ,  m=1,…,M .

The designed scheme is evaluated using 10^6 Monte Carlo runs.

<ref> summarizes the Monte Carlo results for the optimal procedure, along with those of the two-step procedure. The optimal sequential scheme hits the constraints exactly, within the tolerance.
Moreover, the MSPRT achieves smaller empirical detection errors than the constraints, but the estimation constraints are violated since the MSPRT does not take the estimation errors into account.
In <ref>, the empirical run-lengths are summarized for both procedures. Though the two-step procedure has a smaller empirical run-length than the optimal one, this comes at the cost of violating the constraints on the MSE as mentioned previously.





In <ref>, the evolution of both policies over time is shown for three distinct time instances. In this figure, the gray region is the complement of the stopping region of the optimal scheme and the other filled regions are the regions in which the optimal scheme stops and decides for a particular hypothesis.
The hatched areas indicate the regions in which the MSPRT stops and decides in favor of a particular hypothesis.
For n=5 (<ref>), most of the state space corresponds to the complement of the stopping region of the methods and there are four small areas in which both procedures, the optimal and the two-step, stop.
For both methods, these regions are still present for n=10, but larger compared to the ones in <ref>.
In addition, there appear two regions for small/large values of x̅_n in which the optimal procedure stops and decides in favor of _4/_1.
Since a small value of x̅_n increases the certainty about _1, a decision in favor of _4 is not intuitive here.
As we consider a joint detection and estimation problem, the uncertainty about the true hypothesis and the true parameter affect the decision rule. Though a small value of x̅_n implies a high certainty about _1, it leads to a high uncertainty about the true parameter at the same time. Hence, a decision in favor of _4 is less costly.
For the MSPRT, which does not encounter any estimation, this phenomenon is not visible at all.
The four regions for stopping the optimal procedure, which are present at the bottom of <ref> and <ref>, have grown further in <ref>. Though these four regions have grown equally for the optimal method, the regions for stopping the MSPRT and deciding in favor of _1 and _4 now cover a large area of the state space. At the same time, the two regions of the optimal policy, that appeared in <ref>, cover almost the whole state space but make a decision opposite to the MSPRT. 



  

§ CONCLUSIONS

Based on very mild assumptions on the underlying stochastic process, we have developed a flexible framework for sequential joint multiple hypotheses testing and parameter estimation.
The optimal scheme minimizes the expected run-length while fulfilling constraints on the probabilities of falsely rejecting a hypothesis as well as on the MSE
. These procedures have been characterized by a set of non-linear Bellman equations. We have further shown a strong connection of the cost coefficients of the Bellman equations and the performance measures of the method, i.e., the probability of falsely rejecting a hypothesis and the MSE. Based on this connection, we have presented two approaches to obtain the set of optimal cost coefficients. The first approach formulates the problem as a linear program and the second approach uses a projected gradient ascent.
The performance of the optimal procedure has been validated via two numerical examples. The first example has shown the basic properties of the optimal scheme, whereas the second example has been used to show applicability to real world problems. For both examples, Monte Carlo results have been provided. Moreover, the performance gap to a sub-optimal scheme, i.e., a matrix sequential probability ratio test followed by an MMSE estimator, has been shown.
 
sectionappsec










§ PROOF OF COROLLARY <REF>

It has to be shown that if all _m and all _m, m=1,…,M, are finite, then ρ_n+1 is -integrable for all n∈ and all 0≤ n < N.
From the definition of the cost function ρ_n, one can directly see that

    ∫ρ_n+1≤∫ g≤∫ D^⋆_m,n+1 .

With the definition of D^⋆_n+1,m, m=1,…,M, the integral on the right hand side can be written as

    ∫ D^⋆_m,n+1    = ∫_m p(_mn+1) [_m_m, n+1] + ∑_i=1,i≠ m^M _i p(_in+1) 
          = _m ∫ p(_mn, n+1) [_m_m, n, n+1]p(n+1n)n+1
        + ∑_i=1,i≠ m^M _i ∫  p(_in, n+1) p(n+1n)n+1 .

The integral in <ref> can be written as

    p(_mn)∫[_m_m,n,n+1]p(n+1_m,n)n+1 
    
     =   p(_mn)∫∫(_m - [_m_m,n,n+1])^2p(_m,n+1_m,n)n+1_m .

By expanding the square, the double integral in <ref> becomes

    ∫∫_m^2 p(_m,n+1_m,n)n+1_m 
        - 2 ∫∫_m[_m_m,n,n+1] p(_m,n+1_m,n)n+1_m 
        + ∫∫([_m_m,n,n+1])^2 p(_m,n+1_m,n)n+1_m .
    
    =   [ _m^2 _m,n]  - 2 ([_m_m,n])^2 
         + ∫∫([_m_m,n,n+1])^2 p(_m,n+1_m,n)n+1_m .

By applying Jensen's inequality <cit.>, the integral in the last equation can be upper bounded by

    ∫∫([_m_m,n,n+1])^2 p(_m,n+1_m,n)n+1_m 
    ≤   ∫∫[_m^2_m,n,n+1] p(_m,n+1_m,n)n+1_m  
    
     =   [_m^2_m,n] .

Hence, by combining the previous results, it follows for the integral in  <ref> that

    ∫ p(_mnn+1) [_m_m,n, n+1]p(n+1n)n+1≤ 2p(_mn) [_m_m,n] .

Moreover, <ref> reduces to

    ∑_i=1,i≠ m^M _i ∫  p(_in, n+1) p(n+1n)n+1 = ∑_i=1,i≠ m^M _i p(_in) .

Hence, we can conclude that

    ∫ρ_n+1   ≤∫ D^⋆_m,n+1
       ≤_m 2p(_mn)[_m_m, n] + ∑_i=1,i≠ m^M _i p(_in) < ∞ ,

which is finite as long as all _m, _m, m=1,…,M, are finite since the posterior probabilities of _m are finite by definition and the posterior variance is finite by assumption <ref>. 



§ PROOF OF LEMMA <REF>

In order to prove <ref>, we first transfer the cost functions defined in <ref>, the boundary of the stopping region defined in <ref> and the probability measure defined in <ref> to another domain. In the new domain they depend on the sufficient statistic as well as on the posterior probabilities. The posterior probability of _m, m=1,…,M, is denoted by m in what follows. The posterior probabilities are collected in the tuple =(1,…,M), which is defined on the metric state space (, ).
Hence, the combined cost for deciding in favor of _m is given by

    m,n(n,) = ∑_i=1,i≠ m^M _i i + _m m[_m_m, n] .

We can now rewrite the overall cost function as

    ρ̃_n(n,)     = min{g̃(n, ), d̃_n(n, ) }    n<N ,
    ρ̃_N(N,[N])     = g̃(N, [N]) ,

where the cost functions for continuing and stopping the test are defined as

    g̃(n, )    =  min{1,n(n,), …, M,n(n,)} , 
    d̃_n(n, )    =  1 + ∫ρ̃_n+1((), n) p(n) .

The transition kernel of the posterior probabilities is given by [n+1] = n. The equivalent of the probability measure defined in <ref> is given by

    (B×B̃) = P({∈: ()∈ B, _nn∈B̃}) ,

for all elements B of the Borel σ-algebra on  and all elements B̃ of the σ-algebra on . Finally, the counterpart of the boundary of the stopping region defined in <ref> is given by

    = {(n, )∈×: g̃(n, ) = d̃_n(n, ) } .

Before we can prove that the boundary of the stopping region is a P-null set, two auxiliary lemmas have to be introduced.
 
 Let a=(a_1,…,a_M) and let a· denote the element-wise product. Then for all a∈^M, all n∈ and all ∈ it holds that
 
    min{a_1,…,a_M,1}g̃(n,) ≤g̃(n,𝐚·) ≤max{a_1,…,a_M,1}g̃(n,) .




 Since the proof for the upper and lower bound do not differ significantly, only the proof for the lower bound is outlined here. Let a^⋆=min{a_1,…,a_M}, then it holds that
 
    g̃(n, 𝐚·)     = min{1,n(n,𝐚·), …, M,n(n,𝐚·)}
        = a^⋆min{1/a^⋆1,n(n,𝐚·), …, 1/a^⋆M,n(n,𝐚·)} .

It further holds that

    1/a^⋆m,n(n,𝐚·) = ∑_i=1,i≠ m^M a_i/a^⋆_i i + _m a_m/a^⋆m[_m_m, n] ≥m,n(n,) ,

since a^⋆≤ a_m, ∀ m∈{1,…,M}.

 Applying <ref> to <ref>, yields
 
    g̃(n, 𝐚·)     = min{1,n(n,𝐚·), …, M,n(n,𝐚·)}
        = a^⋆min{1/a^⋆1,n(n,), …, 1/a^⋆M,n(n,)}
       ≥ a^⋆min{1,n(n,), …, M,n(n,)}
       ≥min{a^⋆,1}min{1,n(n,), …, M,n(n,)} = min{a_1,…,a_M,1}g̃(n, ) ,

 which is the lower bound stated in <ref>.


 
 Let 𝐚=(a_1,…,a_M) and let 𝐚· denote the element-wise product. Then for all 𝐚∈^M, all n∈ and all ∈ it holds that
 
    min{a_1,…,a_M,1}ρ̃_n(n,) ≤ρ̃_n(n,𝐚·) ≤max{a_1,…,a_M,1}ρ̃_n(n,) .



 Since the proofs for the upper and lower bound do not differ significantly, only the proof for the lower bound is outlined here. First of all, it has to be mentioned that the transition kernel relating the posterior probabilities m and [n+1]m is linear in m, i.e.,
 
    [n+1]m = _nn+1n = _n p(n+1_m,n)/p(n+1n) .

The proof is done via induction. Let a^⋆ = min{a_1,…,a_M} and assume that <ref> holds for some 0<n<N. Then, by applying <ref>, it holds for n-1 that

    ρ̃_n-1(n-1,𝐚·[n-1])     = min{g̃(n-1,𝐚·[n-1]), 1 + ∫ρ̃_n([n-1](n), 𝐚·[n-1]nn-1 p(nn-1) n}
        = min{g̃(n-1,𝐚·[n-1]), 1 + ∫ρ̃_n([n-1](n), 𝐚·[n-1]nn-1 p(nn-1) n} .

By applying <ref> and <ref>, one obtains

    ρ̃_n-1(n-1,𝐚·[n-1])    ≥min{g̃(n-1,𝐚·[n-1]), 1 + a^⋆∫ρ̃_n([n-1](n),[n-1]nn-1) p(nn-1) n}
       ≥min{ a^⋆g̃(n-1, [n-1]), 1 + a^⋆∫ρ̃_n([n-1](n),[n-1]nn-1) p(nn-1) n}
        = min{ a^⋆g̃(n-1, [n-1]), 1 + a^⋆∫ρ̃_n[n-1] } .

We can further state that

    min{ a^⋆g̃(n-1, [n-1]), 1 + a^⋆∫ρ̃_n[n-1] }   ≥min{ a^⋆g̃(n-1, [n-1]), a^⋆ + a^⋆∫ρ̃_n[n-1] }
       ≥ a^⋆min{g̃(n-1, [n-1]), 1 +  ∫ρ̃_n[n-1] } .

With these results, we can conclude that <ref> holds for n-1 if it holds for n. The induction basis is given by n=N, where it holds that

    ρ̃_N(N,𝐚·[N]) = g̃_N(N,𝐚·[N]) ≥ a^⋆g̃_N(N, [N]) = a^⋆ρ_N(N,  [N]) .



Now, with the help of <ref> and <ref>, <ref> can be proven easily by contradiction. Assume, that there exists a non-zero probability () that the test hits the boundary of the stopping region with its next update for some n<N. Since the posterior probabilities [N] are assumed to be continuous random variables a n∈ and an interval [^∙, 𝐚·^∙] with 𝐚=(a_1,…,a_M) and a_m>1 for all m=1,…,M have to exist for which the costs for stopping and continuing are equal. Mathematically, this can be written as

    g̃(n,𝐚·[n]) = 1 + ∫ρ̃_n+1  ∀_n ∈ [^∙, 𝐚·^∙] .

With the previous results we can conclude that

    1 + ∫ρ̃_n+1𝐚·    = 1 + ∫ρ̃_n+1((n+1), 𝐚·n+1n) p(n+1n) n+1
        = 1 + ∫ρ̃_n+1((n+1), 𝐚·n+1n) p(n+1n) n+1
       ≥ 1 + a^⋆∫ρ̃_n+1((n+1), n+1n) p(n+1n) n+1
        > a^⋆ + a^⋆∫ρ̃_n+1((n+1), n+1n) p(n+1n) n+1
        = a^⋆(1 + ∫ρ̃_n+1𝐚·)
     = a^⋆ g(n,) ≥ g(𝐚·) ,

where a^⋆=min{a_1,…,a_M}. The first and last inequality are to <ref> and <ref>, respectively. This contradicts the assumption that there exist a n∈ and an interval [^∙, 𝐚·^∙] in which the costs for stopping and continuing the test are equal.
Due to the fact that  and  only differ in their representation, this contradiction also implies that ()=0,  ∀n∈. This concludes the proof. 



§ PROOF OF LEMMA <REF>

Let

    ρ^'_n,_m(n)     = ∂ρ_n(n)/∂_m = 
          g__m^'(n)    for n∈
    ∂/∂_m∫ρ_n+1   for n∈
denote the derivative of ρ_n with respect to _m, which is defined everywhere on ∖.
Assume for now that the order of differentiation and integration can be interchanged, i.e.,

    ∂/∂_m∫ρ_n+1 = ∫ρ_n+1,_m^' .

In general, the derivative of ρ_n with respect to _m can now be written as

    ρ^'_n,_m(n) = g__m^'(n) + ∫ρ_n+1,_m^' .

On the stopping region , it holds that

    g__m^'(n)     = ∑_i=1,i≠ m^Mni p(_mn) .

With the short-hand notation

    z_n^m = p(n_m)/p(n)

and the property p(nn-1)=p(nn-1), we can write

    g__m^'(n)     = ∑_i=1,i≠ m^Mni p(_mn) = ∑_i=1,i≠ m^Mni p(_m)z_n-1^mp(n_m,n-1)/p(nn-1) .

With the use of <ref>, we can further state that

    ∫ρ_n+1,_m^'     = ∑_i=0,i ≠ m^M∫_{∈n+1i} p(_m)z_n^mp(_m,n)  + ∫_{[n+1]}ρ^'_n+1,_m
        = ∑_i=0,i ≠ m^Mp(_m)z_n^m^m(n+1i) + ∫_{[n+1]}ρ^'_n+1,_m
        = p(_m)z_n^m^m(n+1m̅)  + ∫_{[n+1]}ρ^'_n+1,_m .

Hence, the derivative of ρ_n with respect to _m is

    ρ_n,_m^'(n)     = nm̅ p(_m) z^m_n  + ( p(_m)z_n^m^m(n+1m̅) + ∫_{[n+1]}ρ^'_n+1,_m) ,

which is the expression stated in <ref>.
It is left to show that the order of the integration and the differentiation with respect to _m can be interchanged.
According to the differentiation lemma <cit.> this holds if and only if

 
  * The function ρ_n+1(n+1) has to be -integrable for all 0≤ n < N and all n∈.
 
  * The function ρ_n(n) has to be differentiable for all 0≤ n ≤ N and all n∈.
 
  * A set of functions h_n^m and f_n^m has to exist, which are independent of _m and _m, respectively. It must further hold
 
    |ρ^'_n,_m(m) |   ≤ h_n^m(n) ,  ∀ m∈{1,…,M} , 
    |ρ^'_n,_m(m) |   ≤ f_n^m(n) ,  ∀ m∈{1,…,M} .


Condition 1 is true by <ref>, but conditions 2 and 3 have still to be proven. The proof is only carried out for the derivatives with respect to _m, since it can be proven analogously with the derivatives with respect to _m.
Assume that the differentiation lemma holds for some n≥ 1 and some m∈{1,…,M}. It has now to be shown that the differentiation lemma holds for n-1 as well. On the stopping region, the derivative is given by

    ρ_n-1,_m^'(m) = nm̅ p(_mn) ,

which is well defined and bounded.
On the complement of the stopping region, it holds that

    ∂/∂_mρ_n-1(n)    = ∂/∂_m∫ρ_n[n-1] = 
     p(_m)z_n^m[n-1]^m(nm̅) + ∫_[n]ρ^'_n,_m[n-1] .

Due to the assumption that the differentiation lemma already holds for n, we can show that for n∈ it holds that

    |∂/∂_mρ_n-1(n)|     = |
     p(_m)z_n^m[n-1]^m(nm̅) + ∫_[n]ρ^'_n,_m[n-1] |
       ≤|  p(_m)z_n^m[n-1]^m(nm̅)| + |∫_[n]ρ^'_n,_m[n-1] |
        =  p(_m)z_n^m[n-1]^m(nm̅) + |∫_[n]ρ^'_n,_m[n-1] |
        =  p(_m)z_n^m[n-1]^m(nm̅) + ∫_[n] h_n^m [n-1] .

Hence, the derivative is bounded on  as well. For the induction basis ρ_N it holds that

    ρ_N,_m^' = Nm̅ p(_mN) < ∞ ,

and therefore

    |ρ_N,_m^'| = Nm̅ p(_mN) = h_N^m(n) .

This concludes the proof. 




§ PROOF OF THEOREM <REF>
 
Assuming that the optimal policy as stated in <ref> is used, the detection errors can be written as:

    nm(n)     = 
                                    0    for  n∈nm
    
                                    1    for  n∈nm̅
    [n+1m(n+1)_m,n]    for  n∈
The expected value in <ref> can be rewritten as

    [n+1m(n+1)_m,n]     = ∫n+1m^m 
        = ∫_{∈}∑_i=1,i≠ m^M∈n+1m^m  + ∫_{∈[n+1]}n+1m^m 
        =  ^m(n+1m̅) + ∫_{∈[n+1]}n+1m^m .

For n∈nm, it holds that

    ρ_n,_m^'(n)/p(_m)z_n^m =    0/p(_m)z_n^m = 0 = nm(n)
and for n∈nm̅, it further holds that

    ρ_n,_m^'(n)/p(_m)z_n^m =    p(_m)z_n^m/p(_m)z_n^m = 1 = nm(n) ,
which is exactly the detection error on the stopping region.
It is now left to show that nm=ρ_n,_m^'(n)/p(_m)z_n^m also solves <ref> on the complement of the stopping region. For n∈, it holds that

    ρ_n,_m^'(n)/p(_m)z_n^m    = ^m(n+1m̅)+ ∫_{∈[n+1]}ρ_n+1,_m^'(())/p(_m)z_n+1^m p(_m,n)  , 
    ρ_n,_m^'(n)/p(_m)z_n^m    = ^m(n+1m̅) + ∫_{∈[n+1]}ρ_n+1,_m^'(())/p(_m)z_n^m p(_m,n)/p(n) p(_m,n)  ,
    ρ_n,_m^'(n)/p(_m)z_n^m    = ^m(n+1m̅) + 1/p(_m)z_n^m∫_{∈[n+1]}ρ_n+1,_m^'(()) p(n)  , 
    ρ_n,_m^'(n)     = p(_m)z_n^m^m(n+1m̅) + ∫_{∈[n+1]}ρ_n+1,_m^' ,
which is true by <ref>. Hence, <ref> holds for nm. 







§ PROOF OF THEOREM <REF>

It has to be shown that

    max_≥0,≥0  L_[],[](,) =  max_≥0,≥0 {ρ_0(0) -∑_m=1^M p(_m)(_m + _m)}

attains its maximum for some non-negative and finite values of ^⋆, ^⋆ and that the solutions of <ref> and <ref> coincide.

By applying <ref>, one obtains

    ∂/∂_mL_[],[](,)     = p(_m)(m - ) ,  m=1,…,M ,
    ∂/∂_mL_[],[](,)     = p(_m)(m - ) ,  m=1,…,M .

Since all constraints in <ref> are inequality constraints, ^⋆/^⋆ are solutions of <ref> if ^⋆/^⋆ are positive and the corresponding derivative vanishes, or if ^⋆/^⋆ are zero and the corresponding derivative is negative. The first case, i.e., when the gradient vanishes, holds, if


    m    =   and  m  =  ,   m=1,…,M ,
i.e., the constraints are fulfilled with equality. It now has to be shown that these ^⋆,^⋆ are non-negative and finite. This is only outlined for ^⋆, since it can be shown similarly for ^⋆. We first consider the limit

    lim__m→∞∂/∂_mL_[],[](,)     = p(_m)(0 - )<0 ,

which contradicts the fact that an infinitely large _m is a solution of <ref>, since a negative gradient would only result in a maximum if and only if  ^⋆_m=0. Next, the gradient at _m^⋆=0 needs closer inspection, i.e.,

    ∂/∂_mL_[],[](,)|__m^⋆=0 = p(_m)(m - ) .

At this point, it has to be mentioned again that the cost for rejecting hypothesis _m does not only depend on _m but rather on all _i as well as on _j, j∈{1,…,M}∖ m. Hence, even if ^⋆_m=0, i.e., the constraint on the error probability under _m is not enforced with equality, it can still be satisfied implicitly, as a consequence of the remaining constraints.
We distinguish the two cases in which the resulting detection error is smaller and larger than the error constraint , respectively.
In the former case, the gradient is negative. This results, in combination with _m^⋆=0, in an optimum (complementary slackness).
In the case where the resulting detection error is larger than the constraint, the gradient becomes positive. This contradicts the assumption that _m^⋆ is an optimal solution.
Due to the concavity of L_[],[](,) and the fact that an infinitely large _m results in a negative gradient, a positive gradient for _m=0 implies that there exists a positive and finite _m^⋆ such that the gradient vanishes.
That is, the designed sequential scheme fulfills the requirements and is of minimum run-length by definition.
It can now be easily shown that the optimal objective is the average run-length

    L_[],[](^⋆,^⋆)     =  ρ_0(0) - ∑_m=1^M p(_m)(_m^⋆ + _m^⋆) 
        = [τ] + ∑_m=1^M p(_m)(_m^⋆m + _m^⋆m) - ∑_m=1^M p(_m)(_m^⋆ + _m^⋆)
        = [τ] .

This concludes the proof. 




§ DERIVATION OF THE POSTERIOR PROBABILITIES OF THE HYPOTHESES
 
According to Bayes' theorem, we can calculate the posterior probability as

    p(_mn) = p(n_m)p(_m)/p(n) .
The marginal density is given by

    p(n_m)     = ∫ p(n_m,σ^2)p(σ^2)σ^2 
        =  ∫(2πσ^2)^-n/2exp(-∑_k=1^n(x_k-A_m)^2/2σ^2) exp(-/σ^2)^/Γ()σ^2 .
Using the parametrization given in <ref>, we obtain

    p(n_m)    =   (2π)^-n/2^/()^Γ()/Γ()∫(σ^2,)σ^2 
       =  (2π)^-n/2^/()^Γ()/Γ() .

Hence, the posterior probabilities are given by

    p(_mn) =  K p(_m) (2π)^-n/2^/()^Γ()/Γ() ,

where the normalization constant K is given by

    K=(∑_m=1^M p(_mn))^-1 .




§ DERIVATION OF THE POSTERIOR PREDICTIVE
 
The posterior predictive can be factorized as

    p(n) = ∑_m=1^M p(_m,n)p(_mn) .
The conditional posterior predictive can now by calculated as

    p(_m,n)     = ∫ p(_m,σ^2) p(σ^2_m,n)σ^2
        = ∫()^/Γ() (σ^2)^- - 1exp(-/σ^2) (2πσ^2)^-0.5exp(-(-A_m)^2/2σ^2)σ^2 
        = 1/√(2π)()^/Γ()∫  (σ^2)^- - 1.5exp(- + 0.5(-A_m)^2/σ^2) σ^2 .

Substituting

    =  + 0.5   and    =  + 0.5(-A_m)^2 ,
we can rewrite the conditional posterior predictive as

    p(_m,n)     = 1/√(2π)()^/Γ()
       ×∫  (σ^2)^- - 1exp(-/σ^2) σ^2  .

As the integrand in <ref> is the scaled pdf of the inverse Gamma distribution, we simplify <ref> to

    p(_m,n)     = 1/√(2π)()^/Γ()Γ()/()^ .

Inserting <ref> and <ref> into <ref> and using K as defined in <ref>, the posterior predictive is given by

    p(n) = K (2π)^-n+1/2∑_m=1^M p(_m) Γ()/Γ()^/()^ .


 



§ ACKNOWLEDGEMENTS

The authors would like to thank the anonymous reviewers and the editors for their time and effort.
The work of Dominik Reinhard was supported by the German Research Foundation (DFG) under grant number 390542458.
The work of Michael Fauß was supported by the German Research Foundation (DFG) under grant number 424522268.
 
tandfx
authoryear, open=(,close=)
 
