Stefano Marano1, Ali H Sayed2. 1. DIEM, University of Salerno, via Giovanni Paolo II 132, Fisciano SA, I-84084, Italy. 2. EPFL, School of Engineering, Lausanne, CH-1015, Switzerland.
Abstract
This work focuses on the development of a new family of decision-making algorithms for adaptation and learning, which are specifically tailored to decision problems and are constructed by building up on first principles from decision theory. A key observation is that estimation and decision problems are structurally different and, therefore, algorithms that have proven successful for the former need not perform well when adjusted for the latter. Exploiting classical tools from quickest detection, we propose a tailored version of Page's test, referred to as BLLR (barrier log-likelihood ratio) test, and demonstrate its applicability to real-data from the COVID-19 pandemic in Italy. The results illustrate the ability of the design tool to track the different phases of the outbreak.
This work focuses on the development of a new family of decision-making algorithms for adaptation and learning, which are specifically tailored to decision problems and are constructed by building up on first principles from decision theory. A key observation is that estimation and decision problems are structurally different and, therefore, algorithms that have proven successful for the former need not perform well when adjusted for the latter. Exploiting classical tools from quickest detection, we propose a tailored version of Page's test, referred to as BLLR (barrier log-likelihood ratio) test, and demonstrate its applicability to real-data from the COVID-19 pandemic in Italy. The results illustrate the ability of the design tool to track the different phases of the outbreak.
With the increasing availability of streaming data, modern applications of statistical inference have been faced with new challenges, resulting in novel operative modalities. Traditionally, statistical inference has been focused on estimation and decision problems based on observations coming from a homogeneous/stationary source; the task of the system is to learn from that environment. The challenge of modern systems is to infuse the learning task with new online adaptation capabilities, which are critical in the presence of changes in statistical conditions, drifting environmental characteristics and time-varying data acquisition/storage/transmission modalities. As learning and adaptation are essentially contrasting requirements, it is clear that the data deluge has imposed an important change of perspective.In learning and adaptation contexts, estimation problems were considered first, which led to the adoption of the LMS (least-mean-square) algorithm for the adaptation step, because of its well-known adaptation properties. When dealing with decision problems, it appears natural to maintain the LMS protocol as the basic engine due to its simplicity in order to track drifts in the state of nature. This trend is evident in distributed decision systems based on consensus strategies [1], [2], [3], [4], [5], [6], [7] or diffusion strategies [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20]. Among the latter class, the most successful strategy is the ATC (adapt-then-combine) fusion rule. This strategy is conceived in the form of LMS iterates followed by a convex combination of states [15], [16], [17], [18], [19]. Thus, the state-of-the-art of modern decision systems is built on an estimation-based engine, which may represent a challenge for decision systems operating under the online learning and adaptation paradigm. The fundamental question we pose is if the performance of LMS-based decision systems can be improved by replacing their LMS component with algorithms specifically designed for solving decision tasks. Such modification has the potential of prompting new system architectures and design guidelines.
Contribution and scope
In this paper we design alternative decision-making algorithms that are specifically tailored to decision problems, by building up on first principles from decision theory rather than relying directly on the LMS update. Borrowing ideas from quickest detection, we introduce performance indexes that quantify the tradeoff between learning and adaptation and derive closed-form analytical expressions for the system performance. Using these performance figures, we show that the proposed updating rule outperforms the one based on the LMS iteration. In this contribution, no network aspects are considered and the focus is on the operation of a single agent. In addition, we limit our study to the case of two possible states of nature, say and , which are known to the decision maker. This means that at any time epoch, the observations collected by the agent are independently drawn from one of two distributions, but we do not know which, and this underlying distribution is allowed to change at any time according to arbitrary patterns of the kind .The final part of this article is devoted to a real-word application of the proposed decision procedure. The COVID-19 pandemic due to the new SARS-CoV-2 virus is one of the most serious global crisis of the last decades. One key aspect to counteract this threat is to identify growing/shrinking pandemic phases in a timely manner. Our approach is that such a problem, after identifying the growing/shrinking phases with the statistical hypotheses and , respectively, naturally fits into the framework of binary decisions under the learning and adaptation paradigm. We illustrate that the technique developed in this paper is useful in tracking different phases of the COVID-19 pandemic. As an example, the proposed tool is demonstrated on pandemic data from Italy recorded from the beginning of the pandemic in late February 2020, to the end of July 2021.
Notation
Boldface symbols denote random scalars/vectors and normal fonts their realizations and deterministic quantities. For scalars, the time index (or algorithm iteration number) is enclosed in parentheses. Thus, for instance, denotes the random scalar at time . Conversely, in the case of vectors, the time dependence is indicated by a subscript, as, for example, denotes a random vector evaluated at time . Superscript denotes vector transposition. Statistical expectation, variance, and probability operators are denoted by , , and , respectively. They always are computed under the hypothesis in force , and the pertinent subscript , is usually added to the operator symbol.The remaining part of this article is organized as follows. Section 2 discusses the genesis of LMS in decision contexts. The proposed alternative to LMS is presented in Section 3 and its performance is investigated in Section 5 in terms of the criteria discussed in Section 4. Examples using synthetic data are given in Section 6 while an application to COVID-19 pandemic time-series is discussed in Section 7. Section 8 contains conclusive remarks.
Genesis of the LMS algorithm in decision contexts
Let us start by considering an estimation problem. Let be a zero-mean scalar random variable with variance , and a zero-mean random vector with positive-definite covariance matrix . The quantity is unknown while is observed. The goal is to solve the optimization problem , where is a weight vector and represents a cost function that quantifies the penalty incurred when the unknown is replaced by the linear transformation of the observation. One common choice is the quadratic cost function , in which case the solution is given by , and the linear least-mean-square estimator of given is [21, Th. 8.1, p. 142].A recursive solution to the optimization problem with quadratic cost function is provided by the steepest-descent algorithm: set equal to some initialization vector, and iterate as follows:where the step-size is sufficiently small (less than 2 divided by the largest eigenvalue of matrix ), see [21, Th. 8.2, p. 147]. It can be shown that , which makes it possible to rewrite (1) in terms of the gradient vector . The resulting expression is useful when alternative cost functions are used.What is especially relevant in the adaptive framework is the consideration that the quantities and may not be known beforehand and are expected to vary over time. In these situations, assuming that we have access to streaming data in the form of a sequence of realizations of and , a viable alternative to (1) is obtained if we drop the expectation signs and replace the random variables by their current realizations, yielding the following algorithm: set some initial guess,with a sufficiently small . This stochastic gradient approximation (because the true gradient is replaced by a noisy version thereof) is known as the LMS algorithm, see [21, Th. 10.1, p. 166]. The LMS algorithm learns the data statistics and at the same time is able to track statistical drifts, which are essential characteristics for the design of cognitive intelligent inference systems with learning and adaptation properties.We now move from an estimation to a decision context. Suppose , namely and are scalars, and suppose also for all . By assuming independent and identically distributed (IID) data , formal substitution in (2) gives: ,Note that the right-hand side of (3) is a convex combination: . Iterating (3), we get the output of the LMS algorithm in the form:and we haveFrom (5), we see that the output of the algorithm approximates when the number of iterations is sufficiently large and the step-size is . This property, along with the inherent adaptation ability, motivated the use of (3) in decision problems. Indeed, the algorithm formalized in (3) represents the basic building block for the development of adaptation and learning diffusion algorithms over networks faced with decision problems, which has been addressed in a series of papers [15], [16], [17], [18], [19]. In these works, the term “learning” alludes to a decision about the current state of nature, among an ensemble of possible states characterized by known (not learned online) data distributions. The same meaning to “learning” is given in the present paper.It is useful to mention that the LMS has been advocated in the context of continuous inspection schemes and related control charts. As seen in (4), LMS employs exponentially-scaled weights, which is exactly the idea behind the geometric moving average control charts, see [22, Sec. 2.1.2, p. 28] and [23, Sec. 8.1.2, p. 373]. In these contexts, LMS is known under the name of GMA (Geometric Moving Average) or EWMA (Exponentially Weighted Moving Average). We refer to [22], [23] for details.
Proposed algorithm: Barrier LLR
As discussed in the previous section, in the literature of adaptation and learning, decision problems have been approached by exploiting schemes and protocols initially conceived for estimation problems. Since decision and estimation problems are structurally different in many respects, it makes sense to start anew, with the goal of exploring possible alternatives to the LMS component with better performance for decision tasks. We argue that classical quickest detection tools, if properly adapted, are better suited to solve decision problems, than is the LMS algorithm shown in (3). In the following, we advocate the use of repeated Page’s tests to detect multiple changes in the state of nature. The approach is close in spirit to that adopted in the detection of transients, see, e.g., [24], [25], [26], [27] and references therein.Let us start by considering a standard binary decision problem in which IID data are observed and the following binary hypothesis test must be solved:where are the probability density functions (PDFs) of the data under the two hypotheses and , respectively. These PDFs are assumed to exist and are known to the decision maker. We adopt a non-Bayesian framework, in which the hypotheses represent two different probability spaces. In the complementary setting of the Bayesian formulation, we are given a single probability space with a-priori probabilities assigned to the hypotheses and the goal is to online update these priors by exploiting the data stream. See [28] and the references therein for a recent perspective and [29] for the Bayesian counterpart of the quickest detection procedures exploited in the following.It is well-known that under the most popular optimality criteria, the optimal decision maker for (6) exploits the log-likelihood of the data [30], which isIn view of the IID property of the observations, the optimal decision based on vector iswhere is a suitable threshold, chosen according to the desired optimality criterion [30]. Expression (8) can be regarded as a random walk: ,with step . Note that the observation model in (6) is very general and the linearity of (9) is a consequence of (7) and the IID assumption.The following relationships are well known [the argument “” is suppressed when non-essential]:where and
denote the two Kullback-Leibler (KL) distances (or divergences) between the PDFs and [31]. In (10) we have assumed that these KL distances exist and are strictly positive, which implies that and are distinct over a set of nonzero probability.Assumptions. The following assumptions are used throughout the paper. Under , , the random variables are continuous, with finite first- and second-order moments, and their probability distribution admits a density with respect to the usual Lebesgue measure. In addition, , , and , .In (8), we see that the optimal decision maker compares to a threshold the value of a random walk with positive drift under and negative drift under . This optimal learning scheme is not adaptive, as is easily revealed by the following informal arguments. Suppose that the state of nature is for time steps , and then switches to . Observations are IID under each hypothesis, but their distribution is different under the two hypotheses. Assuming , with high probability the decision statistic takes on very large values because the random walk is drifting to for . For the drift is negative but the random walk “starts” at , implying that the time required to approach the negative values that are typical of hypothesis is very large. A straightforward way to prevent from reaching extreme values is to introduce two barriers , as follows: and for :A more compact expression for the iteration in (11) is: andThe lower and upper barriers limit the range of values that takes on, and hence we can tradeoff adaptation and learning by a careful choice of and . Large values favor the learning (decision) ability of the system, while small values favor its adaptation ability. In the following, the decision procedure based on comparing (11) to a threshold will be referred to as the barrier log-likelihood ratio (BLLR) test: the BLLR decision at any arbitrary time epoch isThe BLLR decision procedure (13), with shown in (11) or (12), represents the proposed alternative to the one based on the LMS engine. For easy reference, the LMS and BLLR procedures are summarized in Algorithms 1
and 2, respectively. Note that the computational complexities of the two procedures are comparable. Essentially, at any iteration, only a function evaluation [to obtain from the measurement ] and a weighted sum (for LMS) or a sum followed by two threshold comparisons (BLLR) are required.
Algorithm 1
LMS.
Algorithm 2
BLLR
LMS.BLLRTypical choices for the threshold appearing in (13) are: , which in the unbounded case of corresponds to the maximum likelihood (ML) decision criterion (also corresponding to the Bayesian criterion when the two hypotheses are equally likely); the mid-point threshold ; or the value of for which the probability of deciding under state of nature takes on a prescribed value (false alarm criterion), as in the Neyman-Pearson formulation [30]. Likewise, when considering the LMS iterate (3), a test similar to (13) is used and the threshold is chosen with the same criteria. For LMS, the mid-point threshold is .Suppose we know that is in force for and is in force for all , with IID data under each hypothesis, and the change epoch is unknown. The celebrated Page’s test for quickest detection of a change in the state of nature is obtained from (11) by setting and letting be the decision threshold: once the decision statistic hits the value , the change in the state of nature is declared and the test stops [22]. This reveals that BLLR test (11) is a generalization of Page’s test. Indeed, the BLLR test in (11) can be seen as an infinite sequence of Page’s tests for successively detecting the changes . To see this, let us assume that is true and set . The BLLR test is equivalent to a Page’s test with threshold for detecting the change , followed by a sign-reversed Page’s test initialized at , driven by negative drifts, with threshold , to detect the successive change , and so forth indefinitely.1
In turn, Page’s test can be regarded as a sequence of Wald’s SPRTs (sequential probability ratio tests) [32], [33], which reveals that the decision algorithm shown in (11) and (13) is a modified version of a sequence of SPRTs. Not surprisingly, the performance analysis of BLLR relies on standard results of sequential analysis, some results of which are collected in Appendix A–Appendix C, for self-consistency.In view of the analogy with sequential analysis, our approach is close in spirit to the SPRT approach pursued in [34] for cooperative sensing. As done in [34], in Section 7 we resort to the GLRT (generalized likelihood ratio test) approach to deal with the presence of unknown parameters. However, the nature of these parameters and the corresponding estimates are structurally different from those in [34], resulting in substantially different decision procedures.Introducing barriers to avoid extreme values of the likelihood also appears in the development of robust approaches to signal detection, as pioneered by Huber in 1965 [35], see, e.g., [30, Sect. III.E.2]. Indeed, the limiting barriers in (11) can be regarded as a soft-limiter nonlinearity and the structure of robust detectors often involves suitably-designed, possibly time-varying, nonlinearities [36]. However, robust detection is concerned with uncertainly in the statistical model of the data, which is a different problem. Note that in our case the nonlinearity is applied to the random walk iterations not to the data, see (11). Besides, time-varying barriers have also been proposed to implement nonparametric/truncated SPRTs [24], [36]; again, these formulations address problems different from that under consideration.
Performance assessment
The performance of quickest detectors are usually expressed in terms of mean time between false alarms and mean decision delay, see, e.g., [22]. In the following, we introduce modified versions of these quality indexes that are better suited to the learning and adaptation scenario under consideration. The analysis exploits standard tools from quickest detection and sequential analysis [22], [23], [32], [33], [37], and results for random walks evolving in a strip [38], [39].
Performance criteria
Performance for BLLR test
We introduce two performance indexes: the error rate
, related to the learning capability, and the expected delay
that quantifies the adaptation capability.Let us consider the learning aspects first. Perhaps, the most natural performance figures would be the probability that exceeds under , and the probability that goes below under . In general, these steady-state probabilities are guaranteed to exist [38], [39], however they are not easy to compute and do not lead to simple closed-form expressions from which insights can be easily gained. We instead introduce performance figures whose computation is tractable. For , consider the following quantity, defined with the state of nature held fixed:2
wherein the sign applies if , and applies if . The quantity in (14) represents the expected time to reach the value starting from , under hypothesis .Using (14), the learning ability of the system is quantified by the two indexesThe interpretation is as follows. The quantity represents the expected time to cross the threshold , yielding a decision in favor of , when the system operates in steady-state and the BLLR iteration is initialized to . The value is referred to as the “typical” value taken by the statistic under . Likewise, represents the expected time to cross the threshold, yielding the decision, when the system operates in steady-state and , where is the “typical” value under . Note that these quantities are related to — but different from — the expected time between false alarms and miss detections, respectively. The expected error time is defined in terms of the quantities in (15), by , and the error index quantifying the learning ability is its inverse, which we call the rate:The second performance index quantifies the adaptation ability and is again defined in terms of shown in (14). Specifically, we consider:For the decision statistic initialized at , represents the expected time needed to hit for the first time the barrier , under steady-state. Likewise, is the expected time for the decision statistic , taking value at epoch 0, to hit for the first time the barrier , under steady-state. The expected delay is defined as the arithmetic meanIn the previous discussion, the “typical” value of the statistic under is , and the “typical” value under is . These choices are natural because and are barriers. With these choices, as we shall see soon, we obtain simple closed-form expressions for the operational characteristic of the decision-maker. However, when comparing the performance of BLLR with that of the LMS test, sensible performance indexes for the BLLR decision-maker are obtained by replacing in (15)–(18) the “typical” values of the decision statistic under the two hypotheses, by the corresponding expected values:wherein we define . The distribution of is investigated, e.g., in [38].
Performance criteria for LMS test
The performance indexes of the LMS test are defined in a way similar to that of BLLR, with the notable difference that, in absence of barriers, one cannot define the typical values of the statistic under the two hypotheses as done before, and we instead rely upon expected values. To elaborate, assuming a steady-state hypothesis , let us introduce the quantity:wherein the sign applies if , and applies if . Recall that is defined in (3) and note the superscript to distinguish quantities related to the LMS test from the corresponding quantities referring to BLLR.As error performance indexes for LMS, we consider the quantities and . The rationale is obvious. For , when the state of nature is and assuming that the iteration starts from , we compute the expected time required to cross the threshold and therefore decide for the opposite hypothesis . Using and , we define the expected error time as the arithmetic mean of these two quantities, and the error rate as the inverse of :Likewise, introducing and , we define the expected delay as
Average run length for Page’s test
The performance of the BLLR test can be computed by borrowing results from the analysis of Page’s test. To show this, it is convenient to introduce a version of the iteration (12) with arbitrary starting point and a single lower barrier at 0. This is exactly the celebrated Page’s test for quickest detection [40]: , andFor , let us define the average run length (ARL):computed under steady-state hypothesis , .In Appendix A, an exact expression for is derived, involving integral equations. Since not much physical insight is gained from these integral representations, we opt for relying on approximate, but simpler and closed-form, performance formulas for . These formulas are derived in Appendix B, exploiting standard Wald’s approximations [32], [33]. The final result is [22, Eq. 5.2.44]:For large , expressions (25) simplify to:For Page’s test, represents the mean time between false alarms and the worst mean delay for detection [22, Eqs. 5.2.18, 5.2.19]. Note the role played by the KL divergences and . The larger and are, the smaller and become, respectively. The former has a positive impact on the performance, the latter has a negative impact. Differently from the classical hypothesis testing problem where an increase of either or both and yields enhanced performance, in quickest detection problem enhanced performance is obtained by increasing and/or , as seen in (25).
Test performance
BLLR test
It is easy to express the performance indexes introduced in Section 4.1.1 in terms of the ARL shown in (25a), (25b). Consider first the case in which the random walk starts at . For the quantities on the left of (15) and (17) we have the obvious equalities:When the random walk starts at , we consider the reversed process obtained by replacing in (11) the sequence of log-likelihoods with its sign-reversed counterpart . Then, for the quantity appearing on the right of (15), we have
where (29a) follows by symmetry; the minus sign “” appended to the ARL in (29b) refers to a “reversed” Page’s test in which the sequence appearing in (23) is replaced by ; and (29c) follows by noting that the ARL is the same of the ARL for a standard Page’s test evolving under with steps whose expectation is . Similar arguments lead toIn the case that the threshold is at the midpoint between the barriers, , the performance figures in (27)–(30) can be expressed in terms of the range
of the detained random walk. Assuming , recalling the definitions of expected error time and expected delay in (16) and (18), we get
where we have defined the effective divergence between the hypotheses asThe inverse of in (31) is the error rate:For , we obtainyielding a simple expression for the operational characteristic
of the BLLR test:The previous equations are insightful. First, the system performance depends on the statistical distributions of the hypotheses only through the effective divergence in (33), which is a combination of the two KL distances and . For , . If one of the two KL distances is much larger than the other, reduces to twice the smallest. The effective divergence appears in the operational characteristic in (36), which quantifies the fundamental trade-off of the decision procedure. Note that is independent of the barriers and , as long as . As it must be, the function is strictly decreasing and convex in . For a fixed , grows with as long as , while it is a decreasing function of for . This behavior is to be interpreted in light of the comments reported at the end of Section 4.2.
LMS test
Simple closed-form approximations for the test performance, similar to those shown in (27)–(30), are not available in the case of the LMS iteration (3). The technical difficulty is that is not a random walk and the stopped martingale approach illustrated in Appendix B does not apply. However, the performance of LMS can be expressed by Fredholm integral equations similar to those in (A.3a), (A.3b).To show this, we follow the approach of [41] with reference to the error figure defined in Section 4.1.2, see (20). The hypothesis in force is , the iteration is initialized to , and the event of crossing the threshold is considered. At the first step of the iteration, two mutually exclusive events can occur. Either causes a threshold crossing, i.e., , an event whose probability we denote by , or . In the latter case the iteration restarts from and the additional expected run length is given by , where the expectation involves the distribution of conditioned to . (These distributions are computed under , even if not explicitly indicated.) Let denote such conditional distribution, which is related to its unconditional counterpart by for , and otherwise. We obtainThe average run length needed for the iteration with initial value to exceed the threshold , can be computed by solving numerically (37). The numerical solution to (37) used in the examples discussed in Section 6 is motivated by the arguments provided in Appendix C. The quantity can be computed similarly to , while and require to consider the reversed random walk process whose steps are . The details are omitted.
Examples
Gaussian shift-in-mean
Consider the following hypotheses with IID observations : for ,It is easily seen that the log-likelihood iswhere . The PDF of is under and under , where denotes a Gaussian distribution with mean and variance . In the present experiment we assume , and , which is also the midpoint threshold because . For the LMS test, different values of the step-size in the range from to 0.3 are considered. In the case of the BLLR test, with little loss of generality, we set the barriers as follows:In this way, we use a single parameter for both LMS and BLLR decision algorithms, with the meaning that smaller values of imply slower adaptation. Special attention is devoted to the slow-adaptation regime .Figure 1 shows the results of computer simulations for both decision schemes. For BLLR we also show the theoretical performance shown in (32) and (34) (theory), as well as the large- expression given in (36) [theory (large )]. For the LMS decision scheme we also show the performance obtained by solving numerically (37) as detailed in Appendix C (LMS numerical). The figure confirms the accuracy of the theoretical formulas for performance prediction.
Fig. 1
Gaussian example of Section 6.1. Operational characteristic (error rate versus expected delay ) for the BLLR and LMS decision procedures. “BLLR simulation” shows the results of computer experiments involving Monte Carlo runs. “BLLR theory” refers to the theoretical formulas (32) and (34), while “BLLR theory (large )” shows the large- approximation (36). The curve in gray labelled as “BLLR theory (correction)” refers to expressions (16) and (18) wherein and are replaced by the expected values as shown in (19), for a fairer comparison with the LMS scheme. “LMS simulation” shows the results of computer experiments involving Monte Carlo runs, while the curve labelled by “LMS numerical” is obtained by solving numerically (37) as detailed in Appendix C.
Gaussian example of Section 6.1. Operational characteristic (error rate versus expected delay ) for the BLLR and LMS decision procedures. “BLLR simulation” shows the results of computer experiments involving Monte Carlo runs. “BLLR theory” refers to the theoretical formulas (32) and (34), while “BLLR theory (large )” shows the large- approximation (36). The curve in gray labelled as “BLLR theory (correction)” refers to expressions (16) and (18) wherein and are replaced by the expected values as shown in (19), for a fairer comparison with the LMS scheme. “LMS simulation” shows the results of computer experiments involving Monte Carlo runs, while the curve labelled by “LMS numerical” is obtained by solving numerically (37) as detailed in Appendix C.The superiority of the BLLR decision algorithm is evident, at least in the regime of small adaptation (large values of ). However, recall from the discussion in Section 4.1.1 that a fair comparison between the two decision schemes requires to modify the performance indexes of the BLLR as indicated in (19). The expectations shown in (19) have been computed numerically and the resulting operational characteristic is shown by the gray curve in Fig. 1, labelled as “theory (correction)”. The substantial superiority of BLLR in the low-adaptation regime is confirmed: for large , we see that the rate scales exponentially with the delay for both the decision schemes, but the exponent is substantially larger for the BLLR decision algorithm.
Example with Gamma distributions
Recall the definition of the Gamma function: , with . Let . With slight abuse of notation we use the symbol to signify that is a Gamma-distributed random variable, whose PDF isFor , it follows that , which is called log-Gamma distribution, having the following PDF:For , we have , where is known as digamma function. We now consider the following hypotheses with IID observations and :Simple algebra shows that the corresponding log-likelihood is , and therefore, with obvious notation:yieldingIn this experiment the barriers for the BLLR decision scheme arewhere is the step-size of the LMS algorithm, and the threshold is . We assume , and .The results of computer experiments for the BLLR and the LMS decision algorithms are shown in Fig. 2
. The comments are similar to those of the Gaussian example. In a nutshell: the formulas for performance prediction are very accurate and the BLLR algorithm outperforms LMS, at least in the small-adaptation regime of large delays.
Fig. 2
Example with Gamma distributions, see Section 6.2. The operational characteristic versus for the BLLR and LMS decision procedures are shown. See the caption to Fig. 1 for details.
Example with Gamma distributions, see Section 6.2. The operational characteristic versus for the BLLR and LMS decision procedures are shown. See the caption to Fig. 1 for details.
Exponentially distributed observations
Our last example involves exponentially distributed observations: with PDF , for and . The two hypotheses arewith . The corresponding log-likelihood is . Defining , for the PDFs of the likelihood one getsand we have: and . In this experiment we assume and . As for the Gamma example, the barriers for the BLLR decision scheme are , and , where is the step-size of the LMS algorithm, and we use the mid-point threshold .The results of computer simulations compared to the theoretical formulas are shown in Fig. 3
. The comments are similar to the previous cases, but in the exponential case the theoretical formulas are less accurate. The slope of the operational curve seems correctly predicted by the analytical formulas, but a multiplicative correction appears to be necessary. This is a manifestation of the poor accuracy of Wald’s approximations of neglecting the excess over the boundaries, which have been exploited to derive the theoretical formulas. Improvements in this regard are possible, e.g., via nonlinear renewal theory [23, Sec. 2.6], [37], but not pursued here. In addition, in the exponential case, the theoretical operational characteristic of the LMS decision scheme is not reported because of instability of the numerical procedure detailed in Appendix C to solve (37).
Fig. 3
Example with exponential distributions, see Section 6.3. The operational characteristic versus for the BLLR and LMS decision procedures are shown. See the caption to Fig. 1 for details.
Example with exponential distributions, see Section 6.3. The operational characteristic versus for the BLLR and LMS decision procedures are shown. See the caption to Fig. 1 for details.
An application related to COVID-19 pandemic
The COVID-19 pandemic, caused by the appearance in late 2019 of a new coronavirus known as SARS-CoV-2, is one of the most serious global crises of the last decades. The multifaceted nature of the threat demands for an extraordinary effort from the scientific community to provide support to informed and rational decision making. In such scenario, signal processing may play an important role and several contributions from this community recently appeared in the literature, including [42], [43], [44], [45] co-authored by one of the present authors; please, see the references therein for useful entry-points to the topical literature. In particular, in [44], [45], a quickest detection procedure — called MAST — is developed. MAST is aimed at detecting the onset of a pandemic phase characterized by an exponential growth of the contagion. The BLLR algorithm developed in the previous sections can be regarded as a repeated application of the MAST to successively detect the growing/shrinking pandemic phases. One difference is that the MAST is based on a surrogate of the log-likelihood (7), to comply with unknown parameters characterizing the statistics of real-world COVID-19 data. In the following, the analogous “GLRT” version of the BLLR is developed and its learning and adaptation properties are demonstrated on COVID-19 time-series and synthetic data.Let us start by considering the classical SIR model of pandemic evolution introduced by Kermack and McKendrick in 1927 [46]. Let3
, and denote the fractions of susceptible, infected, and recovered (or dead) individuals, respectively. Let be the infection rate (infected individuals per unit time), and the recovering rate. The celebrated SIR equations are [47]
with the initial conditions , . We assume , where represents the small fraction of the total population from which the infection originates. To motivate our model, let us focus on the situation in which the fraction of susceptible individuals can be considered approximately constant , implying that the second equation in (49) reduces toSince data about the infections are typically collected on a daily basis, consider a discrete-time version of (50) with unit-step discretization (we loosely use the same notation for the time-discrete version):
for some . From (51) we see that the ratio is constant. It is evident that the real-world phenomenon is governed by “random” versions of the previous deterministic equations. Accordingly, we model the ratio as a random variable. Precisely, we assume:where is a sequence of independent random variables. In light of (53), is referred to as growth rate.By exploiting the publicly available data of the COVID-19 illness spread in Italy (freely downloadable at https://github.com/pcm-dpc/COVID-19/), we obtain4
the sequence of growth rates and verify that the ’s are well-represented by Gaussian random variables; details can be found in [44], [45]. The expected value is time-varying and unknown, and characterizes the specific phase of the pandemic: when the pandemic is under control — a situation here referred to as hypothesis — we have . Conversely, under the alternative hypothesis , and the contagion grows exponentially fast. Upon observing the growth rate sequence , the BLLR machinery can be exploited to online tracking the successive pandemic phases, by deciding between and at any time instant .Lacking knowledge of the sequences of the mean values and , one cannot compute the log-likelihood in (7) and the related BLLR statistic in (11). We then resort to a GLRT approach. As it is well-known, the GLRT approach consists in replacing the unknown parameters appearing in the data distributions with their ML estimates, computed by the data themselves separately for the two hypotheses and the ; details can be found in [30], [48], [49]. The approach is similar to that pursued in [34] in the context of SPRT problems, but the estimates are structurally different. In our case the number of unknown parameters grows linearly in time and the estimates are constrained to and , for all . It is simple to see that
Computing the log-likelihood yields:Substituting the ML estimates (54) and (55) in place of and in (56), one obtainsUsing shown in (57), in place of (7), we get the following “GLRT” version of BLLR: andClearly, the definition of in (57) can also be used in the LMS iteration (3) to obtain a “GLRT” version of LMS. These “GLRT” versions are exactly as described in Algorithms 2 and 1, but take a different sequence in input.We now demonstrate the GLRT versions of BLLR and LMS on COVID-19 time-series recorded in Italy. Top panel of Fig. 4
shows the growth rate sequence , where denotes the index of the day corresponding to the date on the abscissa, from February 25, to July 18, 2021. Bottom panel shows the BLLR and LMS decision statistics computed from . The standard deviation appearing in (57) is set to . For BLLR, we use ; for LMS, . We assume zero threshold for both decision statistics. In the case of BLLR, the first passage is declared on April 13, 2020, followed by passage on July 19, a further on November 27, and so forth. By inspection of Fig. 4, the effectiveness in tracking pandemic phases seems confirmed. For LMS, the first occurrence of is declared on May 4, followed by passage on July 25, followed by on December 20, and so forth. We see that BLLR tends to detect earlier the onset of the pandemic phase, namely reacts more promptly to the change of state of nature. Thus, BLLR appears to be a more sensible decision statistic with respect to LMS, even in terms of GLRT versions for real-world pandemic data. These conclusions are in line with the results of Section 6, and confirmed by the following performance analysis.
Fig. 4
Top: The sequence of growth rates of new positive individuals in Italy, during the COVID-19 pandemic. Bottom: Decision statistics obtained by running the BLLR and LMS procedures on the sequence . Decision is for when the statistic is positive and for otherwise. We use , and .
Top: The sequence of growth rates of new positive individuals in Italy, during the COVID-19 pandemic. Bottom: Decision statistics obtained by running the BLLR and LMS procedures on the sequence . Decision is for when the statistic is positive and for otherwise. We use , and .To investigate the performance of BLLR and LMS on COVID-19 pandemic data, we compute and by computer experiments, assuming that the mean values of the data are IID random variables uniform in under , and uniform in under , for some . Given the mean value, Gaussian data are generated with standard deviation . In the case of BLLR, we explore various values of . In the case of LMS, we compute numerically , and explore various values of . In both cases, Monte Carlo runs are used and the mid-point threshold 0 is selected. The resulting operational characteristics are depicted in Fig. 5
. For both decision schemes, the performance worsens for smaller values of because the growing rate stays closer to unity in both pandemic phases, which makes it more difficult to distinguish between them. We also see that BLLR outperforms LMS, consistently with the results of Section 6, and its superiority is enhanced in the small-adaptation regime of large delays.
Fig. 5
Operational characteristics of the “GLRT” versions of BLLR and LMS, relevant to COVID-19 pandemic control. The mean values of the growth rate are drawn uniformly from under and uniformly from under . We use Monte Carlo runs for each value of (BLLR) and (LMS), mid-point (zero) threshold and .
Operational characteristics of the “GLRT” versions of BLLR and LMS, relevant to COVID-19 pandemic control. The mean values of the growth rate are drawn uniformly from under and uniformly from under . We use Monte Carlo runs for each value of (BLLR) and (LMS), mid-point (zero) threshold and .
Conclusion
We propose a novel algorithm, called BLLR, for learning and adaptation in decision problems. The BLLR decision statistic consists of the data log-likelihood constrained to lie between pre-assigned thresholds and can be interpreted as an infinite sequence of Page’s tests. Performance analysis reveals that BLLR can outperform state-of-the-art LMS-based algorithms, in the slow-adaptation regime. Application to COVID-19 pandemic data is discussed, using a “GLRT” version of BLLR. The proposed approach effectively tracks changes of pandemic phases, providing a rigorous tool to quickly detect the passage from a controlled regime in which the number of new positives tends to decrease or be stable, to a critical regime of pandemic explosion, and vice versa.With the “GLRT” relaxation of BLLR, the statistical distribution of the data may contain unknown parameters. We believe that further relaxations of the formulation adopted in this paper, aimed at addressing non-parametric decision problems, is an interesting direction for future research. Also relevant to big-data applications is a data-centric approach with the decision statistic learned from previously-recorded datasets. Formulations in terms of particle filtering can be useful to address these generalizations. The extension of the proposed method to multi-hypothesis testing also deserves further analysis. Finally, our study is limited to a single decision maker and paves the way for further investigations aimed at designing the diffusion step for a network of interconnected decision makers, much in the same way as the ATC diffusion rule has been advocated to be used in combination with LMS.