Sequential hypothesis testing; changepoint detection; CUSUM; Generalized likelihood ratio; Nonlinear renewal theory; Clinical trial design.

1 INTRODUCTION↩︎

This article provides an overview of Tze Leung Lai’s major contributions to theoretical statistics, biostatistics, applied probability, and sequential analysis. His pioneering work in sequential analysis spans hypothesis testing, changepoint detection, and nonlinear renewal theory. He established fundamental optimality results for the sequential probability ratio test and its extensions, and developed a general framework for testing composite hypotheses. In changepoint detection, he introduced new optimality criteria and computationally efficient procedures that continue to shape the field. He applied many of these and related tools to problems in biostatistics. We review these advances in the broader context of sequential analysis.

To begin, and to offer a sense of Lai’s remarkable life and personality, we present a brief biography.

1.1 Biographical Sketch↩︎

Tze Leung Lai (June 28, 1945 — May 21, 2023) was a pioneering statistician whose influential research profoundly advanced sequential analysis, stochastic modeling, and statistical decision theory. Born in Hong Kong, he earned his B.A. in Mathematics from the University of Hong Kong in 1967 and his Ph.D. in Statistics from Columbia University in 1971 under the supervision of David Siegmund. His doctoral work on sequential estimation and asymptotic optimality laid the foundation for a distinguished career that combined deep theoretical insight with broad practical impact. In recognition of his early and sustained contributions to the field, he received the COPSS Presidents’ Award in 1983 – one of the highest honors in statistics.

Lai made seminal contributions to sequential analysis, extending Wald’s classical framework to adaptive, dependent, and high-dimensional settings. He developed asymptotically optimal stopping rules, sequential confidence procedures, and generalized likelihood ratio (GLR) tests that remain fundamental in modern sequential inference. His unified treatment of sequential testing, estimation, and changepoint detection established enduring theoretical principles that continue to influence online learning, sequential experimentation, and real-time statistical monitoring.

At Stanford University, where he served as the Ray Lyman Wilbur Professor of Statistics, Prof. Lai played a central role in advancing interdisciplinary research. He co-founded and co-directed the Financial and Risk Modeling Institute, where he promoted methodological innovation in financial econometrics, risk management, and financial technology (FinTech). He also directed Stanford’s Center for Innovative Study Design, contributing to the development of adaptive and sequential methods in biostatistics, clinical trials, and health data analysis. His work exemplified the synthesis of rigorous statistical theory with practical applications across science, engineering, and finance.

Throughout his remarkable career, Prof. Lai published more than 300 papers and mentored numerous students and collaborators worldwide. His innovative research reshaped sequential analysis, extending it from its classical foundations into a versatile framework for adaptive learning and data-driven decision-making under uncertainty. He supervised 79 doctoral students and 7 postdoctoral researchers, leaving a lasting legacy of scholarship, mentorship, and inspiration. Figure 1 shows Prof. Lai at two stages of his career.

Prof. Lai married Letitia Chow in 1975. He is survived by Letitia, their two sons, Peter and David, and two grandchildren. He will be remembered for his profound intellect, generosity of spirit, and enduring influence on the field of modern statistics.

1.2 The Remainder of this Article↩︎

The remainder of this article highlights Lai’s influential contributions to sequential hypothesis testing, changepoint detection, and nonlinear renewal theory, as well as selected extensions that continue to inspire new developments.

Section 2 reviews sequential hypothesis testing, covering the near-optimality of the SPRT for two simple hypotheses and Lorden’s 2-SPRT in the modified Kiefer–Weiss problem for general non-i.i.d.settings, as well as the Bayesian and uniform asymptotic optimality of the GLR test with time-varying boundaries for composite hypotheses.

Section 3 surveys Lai’s contributions to changepoint detection and detection–isolation for non-i.i.d.models.

Section 4 discusses Lai’s work on the sequential analysis of multiple change points. Extending Shiryaev’s Bayesian framework, he developed detection rules for unknown pre- and post-change parameters and for multiple change occurrences. Using hidden Markov models with conjugate priors, his methods recursively compute posterior probabilities of change points, enabling efficient detection and estimation.

Section 5 reviews Lai’s foundational results on sequential analysis and nonlinear renewal theory, which provide a unified framework for analyzing boundary-crossing times in stochastic processes—central objects in sequential analysis and related statistical applications. We emphasize asymptotic approximations for stopping times and boundary-crossing probabilities.

Finally, Section 6 reviews some of Lai’s work on biomedical clinical trials, including optimal group sequential designs (Section 6.1) and adaptive designs for mid-trial sample-size re-estimation (Section 6.2).

2 SEQUENTIAL HYPOTHESIS TESTING↩︎

In this section, we review Lai’s contributions to hypothesis testing, particularly his work on establishing the asymptotic optimality of Wald’s SPRT for general non-i.i.d.models, and on the near-optimality of the generalized likelihood ratio test for i.i.d.exponential families in both Bayesian and frequentist settings [1]–[4].

2.1 Near Optimality of Wald’s SPRT for Non-i.i.d.Models↩︎

We begin with the two-hypothesis testing problem for general non-i.i.d.observation models, as addressed by [2]. Let \((\Omega, { \mathscr{F}}, \{{ \mathscr{F}}_n\}, {\mathsf{P}})\), \(n \geqslant 0\), be a filtered probability space, where the sub-\(\sigma\)-algebra \({ \mathscr{F}}_n = \sigma({\mathbf{X}}^n)\) of \({ \mathscr{F}}\) is generated by the sequence of random variables \({\mathbf{X}}^n = \{X_t : 1 \leqslant t \leqslant n\}\) observed up to time \(n\), and \({ \mathscr{F}}_0\) is taken to be trivial. The goal is to test a simple null hypothesis \(\mathop{\mathrm{\mathcal{H}}}_0\!:~ {\mathsf{P}}= {\mathsf{P}}_0\) versus a simple alternative \(\mathop{\mathrm{\mathcal{H}}}_1\!:~ {\mathsf{P}}= {\mathsf{P}}_1\), where \({\mathsf{P}}_0\) and \({\mathsf{P}}_1\) are given probability measures assumed to be locally mutually absolutely continuous. That is, their restrictions \({\mathsf{P}}_0^{n} = {\mathsf{P}}_0|_{{ \mathscr{F}}_n}\) and \({\mathsf{P}}_1^{n} = {\mathsf{P}}_1|_{{ \mathscr{F}}_n}\) to \({ \mathscr{F}}_n\) are equivalent for all \(1 \leqslant n < \infty\). Let \({\mathsf{Q}}^{n}\) denote the restriction of a non-degenerate \(\sigma\)-finite measure \(Q\) on \((\Omega, { \mathscr{F}})\) to \({ \mathscr{F}}_n\).

A sequential test is a pair \(\delta = (T, d)\), where \(T\) is a stopping time with respect to the filtration \(\{{ \mathscr{F}}_n\}_{n \geqslant 1}\), and \(d = d({\mathbf{X}}^T)\) is an \({ \mathscr{F}}_T\)-measurable terminal decision function taking values in the set \(\{0,1\}\). Specifically, \(d = i\) means that hypothesis \(\mathop{\mathrm{\mathcal{H}}}_i\) is accepted upon stopping, i.e., \(\{d = i\} = \{T < \infty,~ \delta~ \text{accepts}~\mathop{\mathrm{\mathcal{H}}}_i\}\). Let \(\alpha_0(\delta) = {\mathsf{P}}_0(d = 1)\) denote the Type I error probability (false positive), and \(\alpha_1(\delta) = {\mathsf{P}}_1(d = 0)\) denote the Type II error probability (false negative) of the test \(\delta\).

For \(i=0,1\), define the class of tests \({\mathbb{C}}(\alpha_0, \alpha_1) = \left\{ \delta : \alpha_i(\delta) \leqslant\alpha_i, ~ i=0,1 \right\}\), that is, the set of procedures whose error probabilities \(\alpha_i(\delta)\) do not exceed the prescribed thresholds \(0 < \alpha_i < 1\).

Consider the general non-i.i.d.model, where the observed random variables \(X_1, X_2, \dots\) may be dependent and non-identically distributed. Under \({\mathsf{P}}_i\), the sample \({\mathbf{X}}^n = (X_1, \dots, X_n)\) has a joint density \(p_i({\mathbf{X}}^n)\) with respect to the dominating measure \({\mathsf{Q}}^n\) for all \(n \geqslant 1\), which can be expressed as \[\label{Jointdensnoniid} p_i({\mathbf{X}}^n) = \prod_{t=1}^n f_{i,t}(X_t \mid {\mathbf{X}}^{t-1}), \quad i = 0, 1,\tag{1}\] where \(f_{i,t}(X_t \mid {\mathbf{X}}^{t-1})\) denotes the conditional density of \(X_t\) given the set of past observations \({\mathbf{X}}^{t-1}\) under hypothesis \(\mathop{\mathrm{\mathcal{H}}}_i\).

For \(n\geqslant 1\), define the log-likelihood ratio (LLR) process between the hypotheses \(\mathop{\mathrm{\mathcal{H}}}_1\) and \(\mathop{\mathrm{\mathcal{H}}}_0\) as \[\begin{align} \lambda_n & = \log \frac{{\mathrm{d}}{\mathsf{P}}_1^{n}}{{\mathrm{d}}{\mathsf{P}}_0^{n}}({\mathbf{X}}^n) = \log \frac{p_{1}({\mathbf{X}}^n)}{p_{0}({\mathbf{X}}^n)} = \sum_{t=1}^n \log \left[\frac{f_{1,t}(X_t \mid {\mathbf{X}}^{t-1})}{f_{0,t}(X_t\mid {\mathbf{X}}^{t-1})}\right] . \end{align}\]

[5], [6] introduced the Sequential Probability Ratio Test (SPRT). Letting \(a_0 < 0\) and \(a_1 > 0\) be thresholds, Wald’s SPRT \(\delta_*(a_0, a_1) = (T_*, d_*)\) is defined by \[\label{SPRT} T_*(a_0, a_1) = \inf\left\{n \geqslant 1 : \lambda_n \notin (a_0, a_1)\right\}, \quad d_*(a_0, a_1) = \begin{cases} 1 & \text{if}~ \lambda_{T_*} \geqslant a_1, \\ 0 & \text{if}~ \lambda_{T_*} \leqslant a_0. \end{cases}\tag{2}\]

For \(i = 0, 1\), let \({\mathsf{E}}_i\) denote expectation under hypothesis \(\mathop{\mathrm{\mathcal{H}}}_i\), i.e., under \({\mathsf{P}}_i\).

In the i.i.d.case, when \(f_{i,t}(X_t \mid {\mathbf{X}}^{t-1}) = f_i(X_t)\) in 1 , Wald’s SPRT exhibits a remarkable optimality property: it minimizes both expected sample sizes, \({\mathsf{E}}_0[T]\) and \({\mathsf{E}}_1[T]\), over the entire class of sequential and non-sequential tests \({\mathbb{C}}(\alpha_0, \alpha_1)\), as established by [7].

However, in the non-i.i.d.case—when the log-likelihood ratio \(\lambda_n\) is no longer a random walk—the Wald–Wolfowitz argument breaks down, and their theorem on SPRT optimality can no longer be applied. This limitation holds even in certain i.i.d.settings, such as invariant sequential tests, which fall outside the scope of their original proof.

This limitation motivated [2] to extend the Wald–Wolfowitz result by proving that the SPRT is also first-order asymptotically optimal as \({\alpha_{\rm max}}= \max(\alpha_0, \alpha_1) \to 0\) for general non-i.i.d.models with dependent and non-identically distributed observations, provided that the normalized log-likelihood ratio \(n^{-1} \lambda_n\) converges \(r\)-quickly to finite limits \(I_i\) under \({\mathsf{P}}_i\), for \(i = 0, 1\).

Given the fundamental importance of this result for both theory and applications, we now present a detailed overview of Lai’s 1981 paper.

First, Lai observed that in the i.i.d.case, the normalized LLR \(n^{-1} \lambda_n\) converges almost surely (a.s.) to \(I_1 > 0\) under \({\mathsf{P}}_1\), and to \(I_0 < 0\) under \({\mathsf{P}}_0\), where \(I_i = {\mathsf{E}}_i[\lambda_1]\), \(i = 0,1\).

This observation naturally led to the idea that, in the general non-i.i.d.case, one should assume the following condition: \[\label{ASconvLLR} \frac{1}{n} \lambda_n \xrightarrow[n \to \infty]{{\mathsf{P}}_i\text{-a.s.}} I_i, \quad i = 0, 1,\tag{3}\] where the limits \(I_0 < 0\) and \(I_1 > 0\) are finite.

This allowed Lai to prove the following weak asymptotic optimality result: if the almost sure convergence condition 3 holds and the SPRT \(\delta_* \in {\mathbb{C}}(\alpha_0, \alpha_1)\) has thresholds \(a_0 \sim \log \alpha_1\), \(a_1 \sim |\log \alpha_0|\) as \({\alpha_{\rm max}}= \max(\alpha_0,\alpha_1) \to 0\), then for \(i=0,1\) and every \(\varepsilon \in (0,1)\), \[\inf_{\delta \in {\mathbb{C}}(\alpha_0, \alpha_1)} {\mathsf{P}}_i\!\left(T > \varepsilon\, T_*(a_0,a_1)\right) \;\to\; 1 \quad \text{as } {\alpha_{\rm max}}\to 0.\] See Theorem 1 in [2].

However, this result is not entirely practical. In applications, one is typically interested in the expectation of the stopping time, or more generally, in its higher moments. In fact, almost sure convergence alone does not even guarantee the finiteness of the expected sample size under general non-i.i.d.models.

To address this, Lai used the concept of \(r\)-quick convergence of the log-likelihood ratio process to establish near-optimality with respect to moments of the stopping time distribution. The definition used in our context is as follows.

For \(i=0,1\), \(\varepsilon > 0\) and \(n \geqslant 1\), define the last exit times \(\tau_{i,\varepsilon} = \sup \left\{ n \geqslant 1 : \left| \lambda_n - n I_i \right| > \varepsilon \right\}\), where \(\sup\{\varnothing\} =0\). For \(r > 0\), the LLR process is said to converge \(r\)-quickly to \(I_i\) under \({\mathsf{P}}_i\) if \[\label{rquickLLR} {\mathsf{E}}_i[\tau_{i,\varepsilon}^r] < \infty \quad \text{for all } \varepsilon > 0, \quad i = 0, 1.\tag{4}\]

The following theorem establishes the first-order asymptotic optimality of the SPRT with respect to \(r\)th moments under the \(r\)-quick convergence conditions 4 for the normalized LLR process \(\lambda_n/n\). This is the outstanding result of Lai’s 1981 paper.

The results of Theorem 1 can be extended to the asymptotically non-stationary case when the LLR process is normalized by a non-linear function \(\psi(n)\) that goes to infinity and increases faster than the logarithmic function, for example, a polynomial function \(\psi(t) = t^k\), \(k > 0\). Specifically, if we require that \[\frac{\lambda_n}{\psi(n)} \xrightarrow[n \to \infty]{{\mathsf{P}}_i\text{-}r\text{-quickly}} I_i, \quad i = 0, 1,\] then in place of the asymptotic expressions 5 the following relations hold: \[\begin{align} \inf_{\delta \in {\mathbb{C}}(\alpha_0, \alpha_1)} {\mathsf{E}}_0[T^r] &\sim \left[ \psi^{-1}\!\biggl(\frac{|\log \alpha_1|}{|I_0|}\biggr) \right]^r \sim {\mathsf{E}}_0[T_*^r], \\ \inf_{\delta \in {\mathbb{C}}(\alpha_0, \alpha_1)} {\mathsf{E}}_1[T^r] &\sim \left[ \psi^{-1}\!\biggl(\frac{|\log \alpha_0|}{I_1}\biggr) \right]^r \sim {\mathsf{E}}_1[T_*^r], \end{align}\] where \(\psi^{-1}(t)\) denotes the inverse function of \(\psi(t)\). The details can be found in [8] and [9].

These results also hold under slightly less restrictive conditions for the LLR process. Specifically, when \(\lambda_n/\psi(n)\) converges to \(I_i\) \(r\)-completely under probability measures \({\mathsf{P}}_i\), that is, if \[\lim_{n\to\infty} \sum_{t=n}^\infty t^{r-1} {\mathsf{P}}_i\Bigl(\bigl\lvert \tfrac{\lambda_t}{\psi(t)} - I_i \bigr\rvert > \varepsilon\Bigr) = 0 \quad \text{for every } \varepsilon > 0, \quad i = 0, 1.\] See [10] for details.

2.2 Near Optimality of Lorden’s 2-SPRT for Non-i.i.d.Models↩︎

[11] advanced the modified Kiefer–Weiss problem for two hypotheses \(\mathop{\mathrm{\mathcal{H}}}_i:{\mathsf{P}}={\mathsf{P}}_i\), \(i=0,1\), by introducing the 2-SPRT—a combination of one-sided SPRTs—and proving its third-order asymptotic optimality as \({\alpha_{\rm max}}\to 0\) under the intermediate measure \({\mathsf{P}}_2\). He later described the structure of optimal tests [12]. Extending this setting to dependent, non-i.i.d.observations, [2] showed that the 2-SPRT is first-order asymptotically optimal whenever the normalized LLRs \(n^{-1} \lambda_n^{(i)} = n^{-1} \log[p_2({\mathbf{X}}^n)/p_i({\mathbf{X}}^n)]\) converge \(r\)-quickly to finite limits \(\eta_i \geqslant 0\) under \({\mathsf{P}}_2\), with \(\max\{\eta_0,\eta_1\}>0\).

Next, letting \(a_0 > 0\) and \(a_1 > 0\) be thresholds, Lorden’s 2‐SPRT \(\delta_2(a_0,a_1) = \bigl(T_2(a_0,a_1),\,d_2(a_0,a_1)\bigr)\) is defined as \[\label{2SPRT} \begin{align} T_2(a_0, a_1) &= \inf\Bigl\{\,n \geqslant 1 : \lambda_n^{(0)} \geqslant a_0 \;\text{or}\; \lambda_n^{(1)} \geqslant a_1\Bigr\}, \quad \inf \{\varnothing\}= \infty, \\ d_2(a_0, a_1) &= \begin{cases} 1, & \text{if } \lambda_{T_2}^{(0)} \geqslant a_0,\\ 0, & \text{if } \lambda_{T_2}^{(1)} \geqslant a_1. \end{cases} \end{align}\tag{6}\]

It is easily shown that the error probabilities of the 2‐SPRT satisfy the inequalities \[{\mathsf{P}}_0\bigl(d_2 = 1\bigr) \;\leqslant\; \exp\left\{-a_0\right\} \,{\mathsf{P}}_2\bigl\{d_2 = 1\bigr\}, \quad {\mathsf{P}}_1\bigl(d_2 = 0\bigr) \;\leqslant\; \exp\left\{-a_1\right\} \,{\mathsf{P}}_2\bigl\{d_2 = 0\bigr\},\] which implies that setting \(a_0 = \log \alpha_0^{-1}\) and \(a_1 = \log \alpha_1^{-1}\) ensures the 2‐SPRT belongs to the class \({\mathbb{C}}(\alpha_0, \alpha_1)\).

The following theorem, which follows from Lai’s Theorem 2 and Corollary 2, is the centerpiece of the modified Kiefer–Weiss non‐i.i.d.theory.

2.3 Near Optimality of GLR Sequential Tests for Composite Hypotheses↩︎

For practical purposes, it is preferable to design tests that minimize \({\mathsf{E}}_\theta[T]\) uniformly over all parameters, rather than only at some intermediate, typically least favorable, point. We focus on sequential tests that are approximately uniformly optimal for small error probabilities, or asymptotically Bayesian when the observation cost is small, in the context of composite hypothesis testing

Let \(X_1,X_2,\dots\) be i.i.d.with density \(f_\theta\), \(\theta \in \Theta = \Theta_0 \,{\textstyle\bigcup}\,\Theta_1 \,{\textstyle\bigcup}\,\Theta_{\mathrm{in}}\). We test \(\mathop{\mathrm{\mathcal{H}}}_0\colon \theta \in \Theta_0\) vs.\(\mathop{\mathrm{\mathcal{H}}}_1\colon \theta \in \Theta_1\), with \(\Theta_{\mathrm{in}}\) an indifference zone. The goal is a sequential test \(\delta=(T,d)\) minimizing \({\mathsf{E}}_\theta[T]\) uniformly over \(\Theta\) within \[\label{classcomp2hyp} {\mathbb{C}}(\alpha_0,\alpha_1) = \bigl\{\delta : \sup_{\theta \in \Theta_0} {\mathsf{P}}_\theta(d = 1) \leqslant\alpha_0,\; \sup_{\theta \in \Theta_1} {\mathsf{P}}_\theta(d = 0) \leqslant\alpha_1 \bigr\}, \quad 0<\alpha_i<1.\tag{8}\]

Since a strictly uniformly optimal test does not exist, we consider asymptotics as the error probabilities vanish, in which case first-order asymptotic optimality in the frequentist setting means \[\label{FOopt} \lim_{{\alpha_{\rm max}}\to 0} \frac{\inf_{\delta\in{\mathbb{C}}(\alpha_0,\alpha_1)} {\mathsf{E}}_\theta[T]}{{\mathsf{E}}_\theta[T]} = 1, \quad \theta \in \Theta.\tag{9}\]

In the Bayesian setting with prior \(\pi(\theta)\), cost \(c\), and loss \(L(\theta)\), the integrated risk is \[\rho_c^\pi(\delta) = \int_{\Theta_0} L(\theta) {\mathsf{P}}_\theta(d=1)\,{\mathrm{d}}\pi(\theta) + \int_{\Theta_1} L(\theta) {\mathsf{P}}_\theta(d=0)\,{\mathrm{d}}\pi(\theta) + c \int_\Theta {\mathsf{E}}_\theta[T]\,{\mathrm{d}}\pi(\theta),\] and tests are asymptotically optimal as \(c\to 0\):

• First‐order: \(\inf_\delta \rho_c^\pi(\delta) = \rho_c^\pi(\delta)(1+o(1))\),

• Second‐order: \(\inf_\delta \rho_c^\pi(\delta) = \rho_c^\pi(\delta)+O(c)\),

• Third‐order: \(\inf_\delta \rho_c^\pi(\delta) = \rho_c^\pi(\delta)+o(c)\).

Consider a single‐parameter exponential family \(\{{\mathsf{P}}_\theta : \theta \in \Theta\}\) with densities \[\label{Expfam} \frac{f_\theta(x)}{f_{\tilde{\theta}}(x)} = \exp\bigl\{(\theta-\tilde{\theta})x - [b(\theta)-b(\tilde{\theta})]\bigr\},\tag{10}\] where \(b(\theta)\) is convex and smooth on \(\widetilde{\Theta} \subset \Theta\). We test \(\mathop{\mathrm{\mathcal{H}}}_0: \underline{\theta}\leqslant\theta \leqslant\theta_0\) vs. \(\mathop{\mathrm{\mathcal{H}}}_1: \theta_1 \leqslant\theta \leqslant\overline{\theta}\), with indifference interval \(\Theta_{\mathrm{in}} = (\theta_0,\theta_1)\) of positive width.

The optimal Bayesian test \(\delta_{\rm opt}=(T_{\rm opt},d_{\rm opt})\) is \[T_{\rm opt} = \inf\{n\geqslant 1 : (S_n,n) \in \mathcal{B}_c\}, \qquad d_{\rm opt}=j \text{ if } (S_n,n)\in \mathcal{B}_c^j,~ j=0,1,\] where \(S_n = \sum_{i=1}^n X_i\) and \(\mathcal{B}_c = \mathcal{B}_c^0 \,{\textstyle\bigcup}\,\mathcal{B}_c^1\) is determined numerically.

[15] derived the test \(\hat{\delta}(\hat{\theta})\), where \(\hat{\theta} = \{\hat{\theta}_n\}\) is the maximum likelihood estimator (MLE) of \(\theta\), as an asymptotic solution, as \(c \to 0\), to the Bayesian problem with the 0–1 loss function: \(L(\theta)=0\) if the decision is correct and \(L(\theta)=1\) otherwise. In this setting, the Bayes integrated risk of a sequential test \(\delta\) is \[\label{AvRisk} \rho_c^\pi(\delta) = \int_{\underline{\theta}}^{\theta_0} {\mathsf{P}}_\theta\bigl(d = 1\bigr)\,{\mathrm{d}}\pi(\theta) + \int_{\theta_1}^{\overline{\theta}} {\mathsf{P}}_\theta\bigl(d = 0\bigr)\,{\mathrm{d}}\pi(\theta) + c \int_{\underline{\theta}}^{\overline{\theta}} {\mathsf{E}}_\theta\bigl[T\bigr]\,{\mathrm{d}}\pi(\theta).\tag{11}\]

Denote by \[\lambda_n(\theta, \theta_i) = \sum_{t=1}^n \log \frac{f_{\theta}(X_t)}{f_{\theta_i}(X_t)} = \bigl[\theta\,S_n - n\,b(\theta)\bigr] \;-\; \bigl[\theta_i\,S_n - n\,b(\theta_i)\bigr]\] the log-likelihood ratio between the points \(\theta\) and \(\theta_i\).

Schwarz’s test prescribes stopping sampling at \(\widehat{T}(\hat{\theta}) \;=\; \min\bigl\{\widehat{T}_0(\hat{\theta}),\,\widehat{T}_1(\hat{\theta})\bigr\},\) where, for each \(i = 0, 1\), \[\label{hatTiSchwarz} \begin{align} \widehat{T}_i(\hat{\theta}) & = \inf\biggl\{\,n \geqslant 1 : \lambda_n\bigl(\hat{\theta}_n, \theta_i\bigr) \;\geqslant\; \lvert \log c\rvert\biggr\} , \quad \inf\varnothing = \infty. \end{align}\tag{12}\] The terminal decision rule \(\hat{d}(\hat{\theta})\) of the test \(\hat{\delta}(\hat{\theta}) \;=\; \bigl(\widehat{T}(\hat{\theta}),\,\hat{d}(\hat{\theta})\bigr)\) accepts \(\mathop{\mathrm{\mathcal{H}}}_0\) if \(\hat{\theta}_{\widehat{T}} < \theta^*,\) where \(\theta^*\) satisfies \(I\bigl(\theta^*, \theta_0\bigr) \;=\; I\bigl(\theta^*, \theta_1\bigr),\) and \[I(\theta, \theta_i) \;=\; {\mathsf{E}}_\theta\bigl[\lambda_1(\theta,\theta_i)\bigr] \;=\; (\theta - \theta_i)\,b'(\theta) \;-\; \bigl[b(\theta) - b(\theta_i)\bigr].\]

Note also that \[\label{SchST} \widehat{T}(\hat{\theta}) \;=\; \inf\Bigl\{\,n \geqslant 1 : n\,\max\bigl[I\bigl(\hat{\theta}_n,\theta_0\bigr),\,I\bigl(\hat{\theta}_n,\theta_1\bigr)\bigr] \;\geqslant\; \lvert \log c\rvert \Bigr\}.\tag{13}\]

A significant advancement in Bayesian theory for testing separated hypotheses about the parameter of the one‐parameter exponential family 10 was made by [16]. [16] introduced a family of GLR tests for separated hypotheses in a single‐parameter exponential family. His tests achieve third‐order asymptotic Bayes optimality, \(\rho_c^\pi(\widehat\delta) = \inf_\delta \rho_c^\pi(\delta) + o(c)\) as \(c \to 0\), by using a slightly reduced threshold and adaptive weight functions correcting for overshoots. These modifications make Lorden’s GLR test nearly optimal compared with the standard Schwartz’s GLR.

[3] improved Bayesian testing of two composite hypotheses by letting the indifference interval \(\Delta = \theta_1 - \theta_0\) shrink as the observation cost \(c \to 0\), resulting in a GLR test with adaptive, time-varying boundaries. This unified approach covers both tests with an indifference zone, \(\mathop{\mathrm{\mathcal{H}}}_0: \theta \leqslant\theta_0\) vs.\(\mathop{\mathrm{\mathcal{H}}}_1: \theta \geqslant\theta_1\), and without, \(\mathop{\mathrm{\mathcal{H}}}_0: \theta < \theta_0\) vs.\(\mathop{\mathrm{\mathcal{H}}}_1: \theta > \theta_0\).

Due to importance of Lai’s theory we now provide details.

In the problem with the indifference zone, Lai proposed replacing the constant threshold \(|\log c|\) in Schwarz’s stopping time 13 with a time‐varying boundary \(g(cn)\): \[\label{LaiST} T^*(\hat{\theta}) \;=\; \inf\Bigl\{\,n \geqslant 1 : n\,\max\bigl[I(\hat{\theta}_n,\theta_0),\,I(\hat{\theta}_n,\theta_1)\bigr] \;\geqslant\; \lvert \log g(cn)\rvert \Bigr\},\tag{14}\] where \(g\colon (0,\infty)\to[0,\infty)\) is any function satisfying conditions detailed below.

If we set \(\theta_0 = \theta_1\) in 14 , the stopping rule simplifies to \[\label{LaiST1} T^*(\hat{\theta}) \;=\; \inf\Bigl\{\,n \geqslant 1 : n\,I\bigl(\hat{\theta}_n,\theta_0\bigr) \;\geqslant\; \lvert \log g(cn)\rvert \Bigr\}.\tag{15}\]

Hence, the stopping rules 14 and 15 provide a unified treatment of both testing problems: separated hypotheses with an indifference zone, and one-sided hypotheses without an indifference zone.

Recall that the integrated Bayesian risk \(\rho_c^\pi(\delta)\) was defined in 11 .

The following theorems summarize Lai’s unified hypothesis‐testing theory. The first addresses asymptotic optimality with an indifference zone, the second without. Let \(J(\theta) = \max\{I(\theta, \theta_0), I(\theta, \theta_1)\}\).

Let \(\{w(t)\}_{t\geqslant 0}\) denote a Wiener process with drift \(\mu\) under \({\mathsf{P}}_\mu\), and define the stopping time \(\tau = \inf\{t>0 : |w(t)| \geqslant h_0(t)\},\) where \(h_0\) is a positive function on \((0,\infty)\) satisfying conditions (2.5) and (2.6) in [3].

In addition to these Bayesian results, Lai also established uniform (over all \(\theta\)) asymptotic optimality of his GLR test within the class \({\mathbb{C}}(\alpha_0,\alpha_1)\) of tests with prescribed error probabilities, thereby solving the frequentist problem 9 . Let \[\alpha_0^* \;=\; \sup_{\theta \leqslant\theta_0} {\mathsf{P}}_{\theta}\bigl(d^* = 1\bigr) \;=\; {\mathsf{P}}_{\theta_0}\bigl(\hat{\theta}_{T^*} > \theta^*\bigr), \quad \alpha_1^* \;=\; \sup_{\theta \geqslant\theta_1} {\mathsf{P}}_{\theta}\bigl(d^* = 0\bigr) \;=\; {\mathsf{P}}_{\theta_1}\bigl(\hat{\theta}_{T^*} \leqslant\theta^*\bigr).\]

Lai’s adaptive GLR \(\delta^*\) uses the MLE \(\hat{\theta}_n\) to adjust the stopping boundary \(g(cn)\), achieving near-optimal performance across \(\theta\) for both indifference-zone and one-sided hypotheses. Unlike Lorden’s random MLE-based boundaries [16], Lai’s are deterministic and handle both cases.

Later, [4] extended these adaptive GLR tests to multiparameter exponential families, addressing one‐sided tests of smooth scalar functions of the vector parameter.

3 QUICKEST CHANGE DETECTION AND ISOLATION↩︎

3.1 Lai’s Changepoint Detection Theory for Non–i.i.d.Data↩︎

In many applications, data may change distribution at an unknown time. A variety of examples were discussed in [9], [14].

In sequential changepoint detection, the goal is to detect distributional changes quickly while controlling false alarms—the quickest change (or disorder) detection problem.

The field originated in quality control with Shewhart’s charts [17], and advanced in the 1950s–70s with the optimal and nearly optimal methods of [18], [19], [20], [21], [22].

For the classical i.i.d.case, let \(X_1, X_2, \dots\)be independent with change point \(\nu\), i.e., \(X_1,\dots,X_{\nu-1}\sim F_0\), \(X_\nu,X_{\nu+1},\dots\sim F_1\), where distributions \(F_0,F_1\) have densities \(f_0,f_1\). A detection rule is a stopping time \(T\) at which a change is declared.

Let \({\mathsf{P}}_\nu, {\mathsf{P}}_\infty\) and \({\mathsf{E}}_\nu, {\mathsf{E}}_\infty\) denote probability measures and expectations when the change point \(\nu\) is fixed (\(1\leqslant\nu < \infty\)) and when \(\nu=\infty\) (i.e., no change ever occurs).

Let \(Z_t \;=\; \log[f_1(X_t)/f_0(X_t)]\) denote the LLR for the \(t\)th observation \(X_t\).

The CUSUM test [19] is defined by \[W_n=\max_{1\leqslant\nu\leqslant n}\sum_{t=\nu}^n Z_t,\qquad T_a=\inf\{n\geqslant 1:W_n\geqslant a\}, \quad \inf \varnothing = \infty,\] where \(a>0\). Here \(W_n\) is the log–GLR statistic and satisfies \[W_n=(W_{n-1}+Z_n)^+,\quad n \geqslant 1, \quad W_0=0.\]

[19] evaluated the CUSUM procedure using the average run length (ARL) to false alarm, \(\mathop{\mathrm{\mathsf{ARL2FA}}}(T) = {\mathsf{E}}_\infty[T]\), and the ARL to detection, \(\mathop{\mathrm{\mathsf{ARL}}}(T) = {\mathsf{E}}_1[T]\), in the case where the change occurs at the very beginning. A more informative measure is the conditional detection delay, \({\mathsf{E}}_\nu[T-\nu \mid T>\nu]\), but under the constraint \(\mathop{\mathrm{\mathsf{ARL2FA}}}(T) \geqslant\gamma\) there exists no procedure that minimizes this quantity uniformly over \(\nu\). Therefore, alternative performance criteria are adopted, most notably Bayesian formulations or minimax approaches that treat \(\nu\) as an unknown parameter.

[20] was the first to address the minimax change detection problem. He considered the class of procedures \({\mathbb{C}}_\gamma = \{T : \mathop{\mathrm{\mathsf{ARL2FA}}}(T) \geqslant\gamma\}\) with \(\gamma \geqslant 1\), and measured detection delay by the worst‐case (double‐supremum) delay \[\label{eq:SADD-Lorden-def} {\mathsf{ESEDD}}(T) = \sup_{0 \leqslant\nu < \infty} \operatornamewithlimits{ess\,sup}_{X_1,\dots,X_{\nu-1}} {\mathsf{E}}_\nu\bigl[(T - \nu + 1)^+ \mid X_1,\dots,X_{\nu-1}\bigr].\tag{16}\] His minimax criterion is \(\inf_{T:\,\mathop{\mathrm{\mathsf{ARL2FA}}}(T)\geqslant\gamma} {\mathsf{ESEDD}}(T).\)

Lorden showed that CUSUM \(T_{a_\gamma}\) with \(a_\gamma=\log\gamma\) is in \({\mathbb{C}}_\gamma\) and \[\inf_{T\in{\mathbb{C}}_\gamma}{\mathsf{ESEDD}}(T)\sim \frac{\log\gamma}{I}\sim{\mathsf{ESEDD}}(T_{a_\gamma}),\quad \gamma\to\infty,\] where \(I={\mathsf{E}}_1[Z_1]\). [23] proved its strict minimax optimality for all \(\gamma\geqslant 1\) with \(\mathop{\mathrm{\mathsf{ARL2FA}}}(T_{a_\gamma})=\gamma\). [24] proposed instead minimizing \[{\mathsf{SEDD}}(T)=\sup_{\nu\geqslant 1}{\mathsf{E}}_\nu[T-\nu+1\mid T\geqslant\nu] \quad \text{ over T \in {\mathbb{C}}_\gamma}.\]

An alternative is the Shiryaev–Roberts procedure \[T^*_A=\inf\{n\geqslant 1:R_n\geqslant A\},\quad R_n=(1+R_{n-1})e^{Z_n},\;R_0=0.\]

As shown in [25], the Shiryaev–Roberts (SR) procedure is asymptotically second-order optimal with respect to Pollak’s \({\mathsf{SEDD}}(T)\) criterion; that is, \[\inf_{T \in {\mathbb{C}}_\gamma} {\mathsf{SEDD}}(T) = {\mathsf{SEDD}}(T^*_{A_\gamma}) + O(1) \quad \text{as } \gamma \to \infty.\]

[24] modified the SR procedure by initializing \(R_0\) from its quasi‐stationary distribution under \(\mathop{\mathrm{\mathcal{H}}}_\infty\). The resulting SRP procedure \(T^P_A\) is third‐order asymptotically minimax optimal: \[\inf_{T\in{\mathbb{C}}_\gamma}{\mathsf{SEDD}}(T)={\mathsf{SEDD}}(T^P_{A_\gamma})+o(1),\quad \gamma\to\infty.\] Similarly, [25] proved third‐order optimality for the SR–\(r\) procedure with \(R_0=r\), where \(r\) is a specially designed fixed point.

[26] extended the Lorden–Moustakides theory to general non-i.i.d.models with dependence, where with pre- and post-change conditional densities \(f_{0,t}(X_t\mid {\mathbf{X}}^{t-1})\) and \(f_{1,t}(X_t\mid {\mathbf{X}}^{t-1})\) the CUSUM statistic is \[W_n = \max_{1 \leqslant\nu \leqslant n} \lambda_n^\nu, \qquad \lambda_n^\nu = \sum_{t=\nu}^n \log \frac{f_{1,t}(X_t\mid {\mathbf{X}}^{t-1})}{f_{0,t}(X_t\mid {\mathbf{X}}^{t-1})}.\] Using change-of-measure and the strong law for LLRs, Lai showed that if \(n^{-1}\lambda_n^\nu \to I > 0\) a.s.and regularity conditions hold, the CUSUM procedure remains first-order asymptotically minimax optimal in the non-i.i.d.case [26].

The following theorem provides necessary details.

The same conclusion, of course, holds for Pollak’s less pessimistic maximal expected detection delay measure, \({\mathsf{SEDD}}(T)\), as well as for the Shiryaev–Roberts procedure.

As mentioned above, the conditional detection delay \({\mathsf{E}}_\nu[T - \nu + 1 \mid T \geqslant\nu]\) cannot be minimized simultaneously for all \(\nu\) under the constraint \(\mathop{\mathrm{\mathsf{ARL2FA}}}(T) \geqslant\gamma\). [26] proposed an alternative criterion, based on the maximal unconditional false alarm probability, \(\sup_{k \geqslant 1} {\mathsf{P}}_\infty(k \leqslant T < k+m),\) and replacing the conditional delay with the unconditional expected detection delay, \({\mathsf{EDD}}_\nu(T) = {\mathsf{E}}_\nu[(T-\nu+1)^+].\)

Lai further proposed setting the time window \(m = m_\alpha\) depending on the false alarm probability constraint \(\alpha\), such that \[\label{malpha} \liminf_{\alpha \to 0} \frac{m_\alpha}{|\log \alpha|} \geqslant\frac{1}{I} \quad \text{and} \quad \lim_{\alpha\to0}\frac{\log m_\alpha}{|\log \alpha|} = 0.\tag{19}\] Then, the corresponding class of detection procedures is defined as \[{\mathbb{C}}_\alpha^{m_\alpha} = {\mathbb{C}}(\alpha) = \left\{T : \sup_{k \geqslant 1} {\mathsf{P}}_\infty(k \leqslant T < k + m_\alpha) \leqslant\alpha \right\}, \quad 0 < \alpha < 1.\]

He then showed that if, instead of the essential supremum condition 17 , the following condition holds: \[\label{Cond1noeessup} \lim_{n \to \infty} \sup_{\nu \geqslant 1} {\mathsf{P}}_\nu \left\{ \max_{1 \leqslant t \leqslant n} \frac{1}{n} \lambda_{\nu+t}^\nu \geqslant I(1 + \varepsilon) \right\} = 0 \quad \text{for all} ~ \varepsilon >0,\tag{20}\] then for any \(T \in {\mathbb{C}}(\alpha)\), as \(\alpha \to 0\), \[\label{CPDLBuniform} {\mathsf{EDD}}_\nu(T) \geqslant\left \{\frac{{\mathsf{P}}_\infty(T \geqslant\nu)}{I} + o(1)\right\} |\log \alpha| \quad \text{uniformly in}~ \nu \geqslant 1.\tag{21}\] See Theorem 2 in [26].

Note that condition 20 holds if \(n^{-1} \lambda_{\nu+n}^\nu \to I\) \({\mathsf{P}}_\nu\)-a.s.

In the non-i.i.d.case, the CUSUM recursion generally fails, making computation costly. [26] proposed a window-limited CUSUM (WLCUSUM) over a sliding window of size \(m_\alpha\): \[W_n^m = \max_{n-m_\alpha \leqslant\nu \leqslant n} \sum_{t=\nu}^{n} \log \frac{f_{1,t}(X_t \mid {\mathbf{X}}^{t-1})}{f_{0,t}(X_t \mid {\mathbf{X}}^{t-1})}, \quad T_a^m = \inf\{ n \geqslant m_\alpha : W_n^m \geqslant a \}.\]

Theorem 4 in [26] shows that \(T_a^m\) is first-order asymptotically optimal as \(\alpha \to 0\), attaining the lower bound 21 when \(a = a_\alpha = \log(2 m_\alpha / \alpha)\) and the left-tail condition 18 holds.

In practical applications, the pre-change conditional densities \(f_{0,n}(X_n \mid {\mathbf{X}}^{n-1})\) are often known, while the post-change densities \(f_{1,n}(X_n \mid {\mathbf{X}}^{n-1})\) are rarely known—they typically depend on unknown parameters \(\theta\). Let \(\{ f_{\theta,n}(X_n \mid {\mathbf{X}}^{n-1}) , \theta \in \Theta \}\) be a parametric family of conditional densities, where the pre-change parameter \(\theta_0 \in \Theta\) is known and the post-change parameter \(\theta \in \Theta\) is unknown. Then we have \[f_{0,n}(X_n \mid {\mathbf{X}}^{n-1}) = f_{{\theta_0,n}}(X_n \mid {\mathbf{X}}^{n-1}), \quad f_{1,n}(X_n \mid {\mathbf{X}}^{n-1}) = f_{\theta,n}(X_n \mid {\mathbf{X}}^{n-1}).\]

Let \({\mathsf{P}}_{\nu,\theta}\) and \({\mathsf{E}}_{\nu,\theta}\) denote probability and expectation when the change occurs at \(\nu\) with post-change parameter \(\theta\), and let \({\mathsf{EDD}}_{\nu,\theta}(T) = {\mathsf{E}}_{\nu,\theta}[(T-\nu+1)^+]\). Define the log-likelihood ratio \[\lambda_n^\nu(\theta) = \sum_{t=\nu}^n \log \frac{f_{\theta,t}(X_t \mid {\mathbf{X}}^{t-1})}{f_{{\theta_0,t}}(X_t \mid {\mathbf{X}}^{t-1})},\] assuming \(n^{-1} \lambda_{\nu+n}^\nu(\theta) \to I_\theta > 0\) a.s.under \({\mathsf{P}}_{\nu,\theta}\).

For a prior \(\pi(\theta)\), the mixture likelihood ratio is \[\overline{\Lambda}_n^\nu = \frac{\int_\Theta \prod_{t=\nu}^n f_{\theta,t}(X_t \mid {\mathbf{X}}^{t-1}) \, {\mathrm{d}}\pi(\theta)}{\prod_{t=\nu}^n f_{{\theta_0,t}}(X_t \mid {\mathbf{X}}^{t-1})}, \quad \bar{\lambda}_n^\nu = \log \overline{\Lambda}_n^\nu.\] The window-limited mixture CUSUM is \[\overline{T}_a = \inf\left\{ n \geqslant m_\alpha : \max_{n-m_\alpha \leqslant\nu \leqslant n} \bar{\lambda}_n^\nu \geqslant a \right\}, \quad \inf \varnothing = \infty.\]

By Doob’s inequality, \(\sup_{k\geqslant 1} {\mathsf{P}}_\infty(k \leqslant\overline{T}_a < k+m_\alpha) \leqslant 2 m_\alpha e^{-a}\), so \(a = \log(2 m_\alpha/\alpha)\) ensures \(\overline{T}_a \in {\mathbb{C}}(\alpha)\).

[26] showed that if \(m_\alpha \to \infty\) and \(\log m_\alpha / |\log \alpha| \to 0\) as \(\alpha \to 0\), and if for each \(\varepsilon>0\) there exist a measurable subset \(\Theta_\varepsilon \subseteq \Theta\) and an integer \(n_\varepsilon\) such that \(\pi(\Theta_\varepsilon) > 0\) and the following condition holds: \[\sup_{n \geqslant n_\varepsilon} \sup_{1 \leqslant\nu \leqslant k} \operatornamewithlimits{ess\,sup}{\mathsf{P}}_{\nu,\theta} \left\{ \frac{1}{n} \inf_{\theta \in \Theta_\varepsilon} \lambda_{k+n}^k(\theta) < I_\theta - \varepsilon \,\middle|\, {\mathbf{X}}^{k-1} \right\} \leqslant\varepsilon,\] then \[{\mathsf{EDD}}_{\nu,\theta}(\overline{T}_{a_\alpha}) \leqslant\frac{{\mathsf{P}}_\infty(\overline{T}_{a_\alpha} \geqslant\nu)}{I_\theta} |\log \alpha| (1 + o(1)) \quad \text{uniformly in }\nu.\] Thus, by the lower bound 21 , the WLCUSUM \(\overline{T}_{a_\alpha}\) is asymptotically optimal.

3.2 Lai’s Changepoint Detection–Isolation Theory for Non–i.i.d.Data↩︎

We consider the quickest changepoint detection problem in the multidecision detection–isolation (or classification–identification) setting with \(N\) possible post-change hypotheses. The goal is to design procedures that asymptotically optimize detection delay while controlling false alarms and misidentifications [9], [14], [27].

Let \(X_1, X_2, \dots\) be observations with change at \(\nu\), and let \(\mathop{\mathrm{\mathcal{H}}}_i\), \(i=1,\dots,N\), denote post-change hypotheses. Denote by \({\mathsf{P}}_\nu^i\) and \({\mathsf{E}}_\nu^i\) the measure and expectation under \(\mathop{\mathrm{\mathcal{H}}}_i\), and by \({\mathsf{P}}_\infty\), \({\mathsf{E}}_\infty\) the no-change scenario.

A sequential change detection–isolation rule \(\delta = (T,d)\) consists of a stopping time \(T\) and a terminal decision \(d \in \{1,\dots,N\}\), where \(d=i\) indicates \(\mathop{\mathrm{\mathcal{H}}}_i\) is accepted at \(T\). In the i.i.d.case, observations are independent with pre-change density \(f_0\) and post-change density \(f_i\) under \(\mathop{\mathrm{\mathcal{H}}}_i\).

In the i.i.d.case, [27] introduced the first minimax theory for change detection and isolation, using the ARL to false alarm or false isolation as risk measures. Observations follow \({\mathsf{P}}_1^i\) if \(\nu=1\) and \(\mathop{\mathrm{\mathcal{H}}}_i\) is true, and follow \({\mathsf{P}}_\infty\) under the nominal regime.

Consider the following sequence of alarm times and final decisions \((T_r, d_r)\): \[T_0 = 0 < T_1 < T_2 < \cdots < T_r < \cdots, \quad\text{and}\quad d_1,\, d_2,\, \ldots,\, d_r,\, \ldots,\] where \(T_r\) is the alarm time of the detection–isolation algorithm applied to the sequence \(X_{T_{r-1}+1}, X_{T_{r-1}+2}, \ldots\).

Let \(X_1, X_2, \dots\) be i.i.d.with pre-change density \(f_0\) and post-change density \(f_i\) under hypothesis \(\mathop{\mathrm{\mathcal{H}}}_i\), \(i=1,\dots,N\). The ARL to the first false alarm of type \(j\) is \({\mathsf{E}}_\infty[\inf\{T_r: d_r=j\}]\), and the ARL to false isolation of type \(j\ne \ell\) under \(\mathop{\mathrm{\mathcal{H}}}_\ell\) is \({\mathsf{E}}_1^\ell[\inf\{T_r: d_r=j\}]\). Define the class \[\label{Nik1} {\mathbb{C}}_\gamma = \Bigl\{ \delta=(T,d) : \min_{0\leqslant\ell \leqslant N} \min_{j\ne \ell} {\mathsf{E}}_1^\ell[\inf\{T_r:d_r=j\}] \geqslant\gamma \Bigr\},\tag{22}\] where \({\mathsf{E}}_1^0={\mathsf{E}}_\infty\). The worst-case expected detection delay is \[{\mathsf{ESEDD}}(\delta) = \max_{1\leqslant i \leqslant N}\sup_{\nu\geqslant 1} \operatornamewithlimits{ess\,sup}{\mathsf{E}}_\nu^i[(T-\nu+1)^+ \mid X_1,\dots,X_{\nu-1}],\] and the minimax goal is \(\min_{\delta\in{\mathbb{C}}_\gamma} {\mathsf{ESEDD}}(\delta)\), usually in the asymptotic regime \(\gamma \to \infty\).

[27] proposed a multi-hypothesis GLR test that is first-order asymptotically optimal but computationally expensive. Later, [28], [29] introduced a minimax formulation using the maximal conditional false isolation probability, \(\sup_{\nu\geqslant 1} {\mathsf{P}}_\nu^\ell(d=j\ne \ell \mid T>\nu),\) with ARL to false alarm \(\mathop{\mathrm{\mathsf{ARL2FA}}}(T) = {\mathsf{E}}_\infty[T]\), leading to the class \[{\mathbb{C}}_{\gamma,\alpha} = \Bigl\{ \delta : \mathop{\mathrm{\mathsf{ARL2FA}}}(T) \geqslant\gamma,\; \max_{i\ne j} \sup_{\nu\geqslant 1} {\mathsf{P}}_\nu^i(d=j \mid T>\nu) \leqslant\alpha \Bigr\}.\] He proposed a matrix CUSUM procedure minimizing \(\max_i {\mathsf{SEDD}}_i(T) = \max_i \sup_\nu {\mathsf{E}}_\nu^i[T-\nu+1 \mid T \geqslant\nu]\) asymptotically as \(\gamma \to \infty, \alpha \to 0\). [30] extended this to per-hypothesis constraints \(\alpha_i\), with an efficient matrix-based procedure minimizing \({\mathsf{SEDD}}_i(\delta)\) for all \(i\) as \(\gamma \to \infty\) and \(\max_i \alpha_i \to 0\).

We now turn to the non-i.i.d.setting and discuss Lai’s important contributions to the theory of change detection and isolation in this more general framework. Assume that the observations are dependent and non-identically distributed, with the pre-change conditional density \(f_{0,t}(X_t \mid {\mathbf{X}}^{t-1})\) for \(t < \nu\), and the post-change conditional density \(f_{i,t}(X_t \mid {\mathbf{X}}^{t-1})\) for \(t \geqslant\nu\), when hypothesis \(\mathop{\mathrm{\mathcal{H}}}_i\) is correct (\(i = 1, \ldots, N\)).

[31] extended [27] by imposing maximal local constraints for all \(\nu \geqslant 1\). For \(i \in \{0,\dots,N\}, j \in \{1,\dots,N\}\), define \[\lambda_n^\nu(i,j) = \sum_{t=\nu}^n \log \frac{f_{i,t}(X_t\mid{\mathbf{X}}^{t-1})}{f_{j,t}(X_t\mid{\mathbf{X}}^{t-1})}.\]

The class controlling local false alarm/isolation is \[\begin{align} {\mathbb{C}}_\alpha(m_\alpha) = \Bigl\{ \delta : &~ \sup_{k\geqslant 1} {\mathsf{P}}_\infty(k \leqslant T < k+m_\alpha) \leqslant\alpha m_\alpha, \\[0.5ex] &~ \max_{1 \leqslant i \leqslant N} \sup_{\nu\geqslant 1} {\mathsf{P}}_\nu^i(\nu \leqslant T < \nu+m_\alpha,\, d \ne i) \leqslant\alpha m_\alpha \Bigr\}. \end{align}\]

The procedure \(\delta_a(m_\alpha) = (T_a, d_a)\) has \[\begin{align} T_a & = \inf \Bigl\{ n : \exists 1 \leqslant i \leqslant N,~ \max_{n-m_\alpha \leqslant k \leqslant n} \lambda_n^k(i,0) - \bigl(\max_{\ell\ne i} \max_{n-m_\alpha \leqslant k \leqslant n} \lambda_n^k(\ell,0)\bigr)^+ \geqslant a \Bigr\}, \\ d_a & = \operatornamewithlimits{arg\,max}_{1 \leqslant i \leqslant N} \max_{T_a-m_\alpha \leqslant k \leqslant T_a} \lambda_{T_a}^k(i,0). \end{align}\]

If \(m_\alpha = O(|\log \alpha|)\), \(a = \log(2N/\alpha)\), and certain regularity conditions hold for the log-likelihood ratios, then the procedure \(\delta_a(m_\alpha)\) is uniformly asymptotically optimal in the class \({\mathbb{C}}_\alpha(m_\alpha)\). For details, see Theorem 7 in [31].

4 MULTIPLE CHANGE-POINT SURVEILLANCE AND ESTIMATION↩︎

In this section, we review Lai’s contributions to Bayesian analysis of multiple change-points.

4.1 An Extended Shiryaev’s Rule for Composite Hypotheses↩︎

Suppose that the observations \(\{ X_n \}_{n\geqslant 1}\) are independent and such that \(X_1, \dots, X_{\nu-1}\) are each distributed according to a common density \(f_0(x)\), while \(X_{\nu}, X_{\nu+1}, \dots\) each follows a common density \(f_1(x) \neq f_0(x)\). Moreover, assume that the change point \(\nu\) has a geometric prior distribution with success probability \(p\), that is, \({\mathsf{P}}(\nu=k) = p(1-p)^{k-1}\) for \(k=1,2,\dots\).

We first consider the case of two simple hypotheses, where \(f_1\) and \(f_0\) are fully known. Assigning a loss of \(c\) for each observation taken after the change point \(\nu\) and a loss of \(1\) for a false alarm before \(\nu\), [21], [32] showed, using optimal stopping theory, that the optimal detection procedure compares the posterior probability of a change with a threshold.

This leads to the Shiryaev procedure, which has the form \[\label{shiryaev46rule46equ3} T_{\rm S}(\eta) = \inf \Bigl\{ n \geqslant 1 : R_{n,p} \geqslant\eta \Bigr\}, \qquad R_{n,p} = \sum_{k=1}^n \prod_{i=k}^n \frac{f_1(X_i)}{(1-p)\, f_0(X_i)}.\tag{23}\]

where \(\eta\) is chosen in such a way that the probability of false alarm (PFA) is exactly equal to \(\alpha\). For large values of the threshold \(\eta\), the average delay to detection of the Shiryaev procedure can be approximated to first order as \(\eta \rightarrow \infty\); see [9]. Moreover, as the geometric prior parameter vanishes (\(p \to 0\)), the Shiryaev statistic converges to the Shiryaev–Roberts statistic, and the corresponding Shiryaev–Roberts rule is asymptotically Bayes risk efficient [24].

Note that the Shiryaev procedure extends naturally to dependent data (see, e.g., [33], [34]).

We now consider the case of two composite hypotheses, that is, \(f_1\) and \(f_0\) are only partially known. Suppose that \(f_1\) and \(f_0\) are not known in advance, but belong to a multivariate exponential family \[\label{dist46exp46family} f_{\theta}(X) = \exp \{ \theta' X - \psi (\theta) \}.\tag{24}\]

Let \(\pi\) be a prior density function \[\label{dist46exp46family46prior} \pi(\theta; a_0, \mu_0) = c(a_0, \boldsymbol{\mu}_0)\, \exp\big\{ a_0 \mu_0' \theta - a_0 \psi(\theta) \big\}, \qquad \theta \in \Theta,\tag{25}\]

in which \(\Theta\) is the parameter space and \[\frac{1}{c(a_0, \boldsymbol{\mu}_0)} = \int_{\Theta} \exp \big\{ a_0 \mu_0' \theta - a_0 \psi(\theta) \big\}\, {\mathrm{d}}\theta, \qquad \boldsymbol{\mu}_0 \in (\nabla \psi)(\Theta),\]

and \(\nabla\) denotes the gradient vector of partial derivatives. The posterior density of \(\theta\) given the observations \(X_1, \dots, X_m\) drawn from \(f_{\theta}\) is \[\label{dist46exp46family46post} \pi \!\left( \theta;\, a_0+m,\, \frac{a_0 \mu_0 + \sum_{i=1}^m X_i }{a_0 + m} \right).\tag{26}\]

Therefore, 25 is a conjugate family of priors and the predictive distribution satisfies \[\int_{\Theta} f_{\theta}(X)\, \pi(\theta; a, \mu)\, {\mathrm{d}}\theta = \frac{c(a, \mu)}{c\!\left(a+1, \tfrac{a \mu + X}{a+1}\right)}.\]

Suppose that the parameter \(\theta\) takes the value \(\theta_0\) for \(t<\nu\) and another value \(\theta_1\) for \(t \geqslant\nu\), and that the change-time \(\nu\) and the pre- and post-change values \(\theta_0\) and \(\theta_1\) are unknown. Following [21], [32], one can use the Bayesian approach that assumes \(\nu\) to be geometric with parameter \(p\) but constrained to be larger than \(n_0\) and that \(\theta_0, \theta_1\) are independent, have the same density function 25 and are also independent of \(\nu\).

Let \(\pi_n = {\mathsf{P}}\{ \nu \leqslant n \mid X_1, \dots, X_n \}\). Whereas \(\pi_n\) is a Markov chain in the case of known parameters \(\theta_0\) and \(\theta_1\), it is no longer Markovian in the present setting with unknown pre- and post-change parameters. Consequently, Shiryaev’s rule, which triggers an alarm once \(\pi_n\) exceeds a threshold, is no longer optimal; see [35], who suggests applying dynamic programming to determine the optimal stopping rule, but also notes that “it is generally difficult to determine the optimal stopping boundaries.” Due to this complexity, [36] introduced a more tractable myopic (two-step-ahead) policy in the univariate Bernoulli case (\(X_i \in \{0,1\}\)).

[37] introduced a modification of Shiryaev’s rule and showed it is asymptotically Bayes as \(p\rightarrow 0\). Let \({\cal F}_t\) be the \(\sigma\)-field generated by \(X_1, \dots, X_t\) and \[\pi_{0,0}=c(a_0, \mu_0), \quad \pi_{i,j} = c \Big(a_0 + j-i+1, \frac{a_0 \mu_0 + \sum_{t=i}^j X_t}{a_0+j-i+1} \Big).\]

Note that for \(n_0 < i \leqslant n\), \[\label{exShi46equ1} {\mathsf{P}}\{ \nu = i |{\cal F}_n\} \propto p(1-p)^{i-1} \pi_{0,0}^2 \big/ \pi_{1,i-1} \pi_{i,n}, \quad {\mathsf{P}}\{ \nu >n |{\cal F}_n\} \propto (1-p)^n \pi_{0,0} \big/ \pi_{1,n}.\tag{27}\]

The normalizing constant is determined by the fact that all the probabilities in 27 sum to 1. Let \(p_{i,n}={\mathsf{P}}\{ \nu=i| {\cal F}_n \}\) be the posterior probability given the observed samples up to time \(n\), we then have \[{\mathsf{P}}(n_0 < \nu \leqslant n | {\cal F}_n) = \sum_{i=n_0+1}^n p_{i,n} = \frac{ \sum_{i=n_0+1}^n {\mathsf{P}}(\nu=i | {\cal F}_n) }{ \sum_{i=n_0+1}^n {\mathsf{P}}(\nu=i| {\cal F}_n) + {\mathsf{P}}(\nu>n | {\cal F}_n) },\]

in which \({\mathsf{P}}\{ \nu = i |{\cal F}_n\}\) and \({\mathsf{P}}\{ \nu >n |{\cal F}_n\}\) are given by 27 . Therefore, Shiryaev’s stopping rule in the present setting of unknown pre- and post- change parameters can again be written in the form of 23 with \[R_{n,p}=\sum_{i=n_0+1}^n \frac{\pi_{0,0} \pi_{1,n}}{ (1-p)^{n-i} \pi_{1,i}\pi_{i,n} }.\]

[37] proposed an extended Shiryaev’s rule \[\label{exShiryaev46rule1} T_{\rm exShi}=\inf \{ n>n_p: {\mathsf{P}}(\nu\leqslant n| \nu\geqslant n-k_p, {\cal F}_n) \geqslant\eta_p \},\tag{28}\]

and showed that it is asymptotically optimal as \(p\rightarrow 0\), for suitably chosen \(k_p\), \(\eta_p\) and \(n_p \geqslant n_0\); see [37]. Since \[{\mathsf{P}}(\nu\leqslant n | \nu\geqslant n-k_p, {\cal F}_n) = \frac{ \sum_{i=n-k_p}^n {\mathsf{P}}(\nu=i | {\cal F}_n) }{ \sum_{i=n-k_p}^n {\mathsf{P}}(\nu=i| {\cal F}_n) + {\mathsf{P}}(\nu>n | {\cal F}_n)},\]

we can use 27 to rewrite 28 in the form \[\label{exShiryaev46rule2} T_{\rm exShi}=\inf \Big\{ n > n_p: \sum_{i=n-k_p}^n \frac{\pi_{0,0} \pi_{1,n}}{(1-p)^{n-i+1} \pi_{1,i-1}\pi_{i,n}} \geqslant\gamma_p \Big\}.\tag{29}\]

This has essentially the same form as Shiryaev’s rule 23 with the obvious changes to accommodate the unknown \(f_{\theta_0}\) and \(f_{\theta_1}\), and with the important sliding window modification \(\sum_{i=n-k_p}^n\) of Shiryaev’s sum \(\sum_{i=n_0+1}^n\), which has too many summands to trigger false alarms when \(\theta_0\) and \(\theta_1\) are estimated sequentially from the observations.

4.2 Surveillance of Multiple Change-Points↩︎

Quickest detection and sequential surveillance are closely related but differ in scope and objectives. Quickest detection targets a single change in a process and is widely applied in areas such as fault detection, navigation integrity monitoring, and radar/sonar signal processing. Sequential surveillance, by contrast, addresses multiple changes over time—commonly called multiple change-point detection—with applications in finance, cybersecurity, and public health. A common strategy is detection–segmentation, where each detected change-point becomes the starting point for the next search. However, false alarms or detection delays at one stage can propagate and degrade the performance of subsequent detections.

[38] introduced a Bayesian model for multiple change-points in multivariate exponential family distributions with unknown pre- and post-change parameters—a key setting in sequential surveillance. Assume observations \(X_1,\dots,X_t,\dots\) follow the multivariate exponential family 24 where \(\theta_t\) may change occasionally, with indicators \(I_t = 1_{{\theta_t \neq \theta_{t-1}}}\) that are i.i.d. Bernoulli\((p)\). When \(I_t=1\), the new \(\theta_t\) is sampled from \(\pi\). The conjugacy of the prior yields explicit formulas for both sequential (filtering) estimates \({\mathsf{E}}(\mu_t|{\boldsymbol{X}}^t)\) and fixed-sample (smoothing) estimates \({\mathsf{E}}(\mu_t|{\boldsymbol{X}}^n)\), where \(\mu_t=\nabla\psi(\theta_t)\).

For the above model, [39], [40] proposed a general hidden Markov filtering approach that yields recursive and tractable estimators of multiple change-points and is computationally very efficient. Let \(K_t = \max \{ s\leqslant t: I_s =1 \}\) be the most recent change-time \(K_t\) up to \(t\). Denote by \(p_{it} = P(K_t = i |{\cal Y}_t)\) the probability that the most recent change-time up to time \(t\) is \(K_t\) and \(f(\cdot | \cdot)\) the conditional density. Since \(K_t\) can take values from 1 to \(t\), [40] showed that the posterior distribution of \(\theta_t\) given \({\boldsymbol{X}}^t\) is a mixture of distributions, \[\label{mcp46equ6462} f(\theta_t |{\boldsymbol{X}}^t) = \sum_{i=1}^t p_{it} \pi(\theta_t; a_0+t-i+1, \bar{X}_{i,t} ).\tag{30}\]

where \(\bar{X}_{i,j} = (a_0 \mu_0 + \sum_{k=i}^j X_k) \big/ (a_0 + j-i+1)\), \(j\geqslant i\), is the posterior mean, and the mixture weights \(p_{it}\) can be represented recursively by \[\label{mcp46equ7} p_{it} = \frac{p_{it}^*}{\sum_{j=1}^t p_{jt}^*}, \qquad p_{it}^* = \left\{ \begin{array}{ll} p \pi_{0,0} / \pi_{t,t} & ifi=t, \\ (1-p) p_{i,t-1} \pi_{i,t-1} / \pi_{i,t} & ifi<t, \end{array} \right.\tag{31}\]

in which \(\pi_{0,0}=c(a_0, \mu_0)\) and \(\pi_{i,j} = c(a_0 + j-i+1, \bar{X}_{i,j} )\). Then the change-point probability, the probability that the most recent change-point occurs at one of the times \(\{s, s+1, \dots, t\}\), and the posterior mean at \(t\) are given, respectively, by \[\label{mcp46filter46est} {\mathsf{P}}(I_t=1|{\boldsymbol{X}}^t) = p_{tt}, \quad {\mathsf{P}}( s\leqslant K_t \leqslant t) = \sum_{i=s}^t p_{it}, \quad {\mathsf{E}}(\mu_t | {\boldsymbol{X}}^t) = \sum_{i=1}^t p_{it} \bar{X}_{i,j}.\tag{32}\]

To obtain a surveillance rule, [38] use sliding windows \(\sum_{i=n-k_p}^n\) as in 29 but with the summands modified to be the posterior probabilities \(p_{in}\) that the most recent change-point up to time \(n\) occurs at \(i\). Then the surveillance rule based on the assumption of multiple change-points (MCP) can be expressed as \[\label{mcp46rule2} T_{\rm MCP}= \inf \big\{ n> n_p: \sum_{i=n-k_p}^n p_{in} \geqslant\gamma_p \big\}.\tag{33}\]

[38] showed that, by suitably choosing \(k_p\) and \(\gamma_p\), the surveillance rule 33 achieves asymptotically optimal Bayes and frequentist performance, using the definitions of false alarm rate and detection delay for multiple change-points.

In practice, two main issues arise when implementing the surveillance rule 33 as \(n\) grows: estimating the hyperparameters \(p\), \(a_0\), and \(\mu_0\), and handling computational cost. Since the forward filter \(\theta_t|{\boldsymbol{X}}^t\) depends on these hyperparameters, [40] proposed an empirical Bayes approach. From \(p_{it}\), the likelihood of \((p,a_0,\mu_0)\) is \[\label{mcp46llh} \prod_{t=1}^n f(X_t | {\boldsymbol{X}}_{t-1}) = \prod_{t=1}^n \Big( \sum_{i=1}^t p_{it}^* \Big),\tag{34}\]

Because \(X_t\) are exchangeable with mean \(\mu_0\), we estimate \(\mu_0\) by the sample mean \(\widehat{\mu}=n^{-1}\sum_{t=1}^n X_t\). A simple choice \(a_0=1\) treats \(\widehat{\mu}\) as one pseudo-observation at a change-point. The key parameter is \(p\), the relative frequency of change-points. Substituting \(a_0=1\) and \(\widehat{\mu}\) into 34 , \(p\) can be estimated by maximizing \(l(p)=\sum_{t=1}^n \log (\sum_{i=1}^t p_{it}^*)\), using a grid search over \({2^j/n: j_0 \leqslant j \leqslant j_1}\).

To reduce the linear computational complexity in 31 , [40] proposed the bounded complexity mixture (BCMIX) approximation. The idea is to retain the most recent \(m\) weights, discard the smallest among the remaining ones, and then reweight at each step. Let \({\cal K}_{t-1}\) be the set of indices kept at stage \(t-1\) (with \({\cal K}_{t-1} \supset {t-1, \dots,t-m}\)). At stage \(t\), compute \(p_{i,t}\) as in 31 for \(i \in {t} \,{\textstyle\bigcup}\,{\cal K}_{t-1}\), identify the smallest \(p_{i,t}\) among indices \(\leqslant t-m\), remove it, and define \[{\cal K}_t = \{ t \} \,{\textstyle\bigcup}\,({\cal K}_{t-1} - \{ i_t \}), \quad p_{i,t} = \Big( p_{i,t}^* \Big/ \sum_{j \in {\cal K}_t } p_{j,t}^* \Big), \quad i \in {\cal K}_t.\]

Thus, \(|{\cal K}_t| \leqslant M\), ensuring bounded complexity regardless of \(n\).

The surveillance rule 33 with BCMIX can monitor complex systems undergoing multiple changes, with applications in finance and beyond. For example, [41] applied it to detect shifts in firms’ credit rating migration generators.

4.3 Inference on Multiple Change Points↩︎

The hidden Markov filtering approach underlying the surveillance rule 33 can also be extended to estimate \(\theta_t\) for each \(t=1,\dots,n\) from observations \(X_1,\dots,X_n\). This yields the smoothing estimate of \(\theta_t\) given \({\boldsymbol{X}}^t\). Because the number and locations of change-points are unknown, direct estimation is computationally demanding. By combining the forward and backward filters, however, smoothing estimates can be obtained efficiently.

In particular, [40] showed how to derive the posterior distribution of \(\theta_t \mid {\boldsymbol{X}}^n\) by applying Bayes’ theorem to combine the forward filter \(\theta_t \mid {\boldsymbol{X}}^t\) with the backward filter \(\theta_t \mid {\boldsymbol{X}}^{t+1,n}\). The backward filter is obtained by reversing time and is expressed as \[\label{mcp46equ8} f(\theta_t | {\cal Y}_{t+1,n}) = p \pi(\theta_t; a_0, \mu_0) + (1-p) \sum_{j=t+1}^n q_{j,t+1} \pi (\theta_t ; a_0 + j-t, \bar{X}_{t+1,j}),\tag{35}\]

where \(q_{jt} = q_{jt}^*/ \big( \sum_{l=t}^n q_{lt}^*)\) and \[\label{mcp46equ9} q_{j,t}^* =\left\{ \begin{array}{ll} p \pi_{0,0} / \pi_{t,t} & ifj=t,\\ (1-p) q_{j,t+1} \pi_{t+1,j}/ \pi_{t, j} & ifj>t. \end{array} \right.\tag{36}\]

By Bayes’ theorem, the smoother is proportional to the product of forward and backward filters divided by the prior distribution, that is, \[\label{mcp46bayes} f(\theta_t | {\boldsymbol{X}}^n) \propto f(\theta_t | {\boldsymbol{X}}^t) f(\theta_t | {\boldsymbol{X}}^{t+1,n}) \big/ \pi(\theta; a_0, \mu_0).\tag{37}\]

Combining 35 with 30 , and noting that \[\pi \big(\theta; a_0+t-i+1, \bar{X}_{i,t} \big) \frac{\pi \big(\theta; a_0+j-t, \bar{X}_{t+1,j} \big) }{ \pi \big(\theta; a_0, \mu_0 \big) } = \frac{\pi_{it} \pi_{t+1,j}}{\pi_{ij} \pi_{00} } \pi \big(\theta; a_0+j-i+1, \bar{X}_{ij} \big),\]

one can use 37 to obtain the smoother \(\theta_t | {\boldsymbol{X}}^n\), which is expressed as \[\label{mcp46sm1} f(\theta_t | {\boldsymbol{X}}^n) = \sum_{1\leqslant i \leqslant t \leqslant j \leqslant n} \beta_{ijt} \pi(\theta_t; a_0 + j-i+1, \bar{X}_{i,j}),\tag{38}\]

where \(\beta_{ijt}= \beta_{ijt}^* \big/ P^*_t\), \(P^*_t = p + \sum_{1 \leqslant i \leqslant t < j \leqslant n} \beta_{ijt}^*\), and \[\label{mcp46sm2} \beta_{ijt}^* = \left\{ \begin{array}{ll} p p_{it} & ifi\leqslant t=j, \\ (1-p) p_{it} q_{j,t+1} \pi_{it} \pi_{t+1,j} \big/ \pi_{ij} \pi_{00} & ifi\leqslant t <j. \end{array} \right.\tag{39}\]

From the above, it follows that the change-point probability and posterior mean at time \(t\) are expressed as \[\label{mcp46sm3} {\mathsf{P}}(I_{t+1}=1|{\boldsymbol{X}}^n) = p \big/ P^*_t, \qquad {\mathsf{E}}(\mu_t | {\boldsymbol{X}}^n) = \sum_{1\leqslant i \leqslant t \leqslant j \leqslant n} \beta_{ijt} \bar{X}_{i,j}.\tag{40}\]

Computing the smoothing estimates using 40 for all \(t=1, \dots, n\) results in a computational complexity \(O(n^3)\). [40] showed how to apply the BCMIX idea to the smoother 38 . We first define the forward filter \(\theta_t \mid {\boldsymbol{X}}^t\) as in the preceding section. For the backward filter \(\theta_t \mid {\boldsymbol{X}}^{t+1,n}\), let \(\widetilde{\cal K}_{t+1}\) denote the set of indices retained at stage \(t+1\) (with \(\widetilde{\cal K}_{t+1} \supset {t+1, \dots, t+m}\)). At stage \(t\), compute \(q_{j,t}\) as in 36 for \(j \in {t} \,{\textstyle\bigcup}\,\widetilde{\cal K}{t+1}\), identify the smallest \(q{j,t}\) among indices \(\geqslant t+m\), remove it, and define \[\widetilde{\cal K}_t= \{ t \} \,{\textstyle\bigcup}\,(\widetilde{\cal K}_{t}- \{ j_t \}), \qquad q_{j,t} = \Big( q_{j,t}^* \Big/ \sum_{j \in \widetilde{\cal K}_t } q_{j,t}^* \Big), \quad j \in \widetilde{\cal K}_t.\]

This gives a BCMIX approximation to the backward filter \(\theta_t \mid {\boldsymbol{X}}^{t+1,n}\). Then, for each \(t=1,\dots,n-1\), the BCMIX approximation to the smoother 38 can be obtained by combining the forward and backward BCMIX filters with selected weight indices \({\cal K}_t\) and \(\widetilde{\cal K}_{t+1}\) via Bayes’ theorem: \[f(\theta_t | {\boldsymbol{X}}^n) \approx \sum_{i\in {\cal K}_t, \;j\in \widetilde{\cal K}_{t+1}} \widetilde{\beta}_{ijt} \pi(\theta_t; a_0 +j-i +1, \bar{X}_{i,j}),\]

in which \[\widetilde{\beta}_{ijt}=\beta_{ijt}^* / \widetilde{P}^*_t, \qquad \widetilde{P}^*_t=p + \sum_{1\leqslant t \leqslant n, i\in {\cal K}_t, j \in \widetilde{\cal K}_{t+1} }\beta_{ijt}^*,\]

and \(\beta_{ijt}^*\) is given by 39 for \(i\in {\cal K}_t\) and \(j \in \widetilde{\cal K}_{t+1}\). The BCMIX approximation to the change-point probability and posterior mean at time \(t\) are therefore \[\widehat{{\mathsf{P}}}(I_{t+1}=1 | {\boldsymbol{X}}^n) = \frac{p}{P_t^*}, \qquad \widehat{{\mathsf{E}}}( \mu_t | {\boldsymbol{X}}^n) = \sum_{i\in {\cal K}_t, \;j \in \widetilde{\cal K}_{t+1} } \widetilde{\beta}_{ijt} \bar{X}_{i,j}.\]

The BCMIX approximation greatly reduces the computational cost of smoothing estimates for \(\theta_t\) (\(t=1,\dots,n\)). In multiple changepoint problems with large \(n\) (e.g., \(n \sim 10^6\)), the full sequence \(\{\theta_t\}_{1\leqslant t \leqslant n}\) can be estimated with \(O(n)\) complexity. This efficiency demonstrates that the Markov filtering framework, combined with BCMIX, can be effectively applied to scalable change-point estimation in complex systems.

Hidden Markov models for multiple change points, along with their smoothing estimates, provide a powerful framework for analyzing complex dynamical systems with abrupt shifts. Building on this approach, [42] developed changepoint autoregressive GARCH models to infer discrete- time volatility dynamics and simultaneous changes in autoregressive coefficients and volatilities from asset prices. Similarly, [43] modeled firms’ credit rating transitions as piecewise homogeneous Markov chains with unobserved structural breaks. Beyond finance, these models have been applied to genomic sequence data, including array-based comparative genomic hybridization [44], parent-specific DNA copy number in tumors [45], and chromatin immunoprecipitation sequencing [46].

5 NONLINEAR RENEWAL THEORY↩︎

In this section, we summarize Lai’s contributions to nonlinear renewal theory and to multivariate Markov renewal theory.

5.1 Nonlinear and Markov Nonlinear Renewal Theories↩︎

Motivated by the study of boundary crossing times and their importance in sequential analysis, nonlinear renewal theory was developed in the seminal works [47]–[50]. Comprehensive treatments of the classical approaches are given in the monographs [14], [51], [52]. Building on these foundations, [53] established several general results that broadened the scope of the theory. More recently, the framework has been extended to the multivariate setting by [54], thereby opening new directions for applications.

Let \(X_1, X_2, \ldots\) be i.i.d.random variables with common distribution \(F\) and finite, positive mean \(\mu = {\mathsf{E}}[ X_1], 0 < \mu < \infty;\) and \(S_n = \sum_{k=1}^n X_k,~n \geqslant 1,\) denote the partial sums. Let \(\{Z_n = S_n + \eta_n, n\geqslant 1\}\) be a perturbed random walk in the following sense: \(S_n\) is a random walk, \(\eta_n\) is \({\cal F}_n\)-measurable, where \({\cal F}_n\) is the \(\sigma\)-algebra generated by \(\{S_k, 1 \leqslant k \leqslant n\}\). Let \(\eta_n\) be slowly changing, i.e. \(\frac{1}{n} \max_{1 \leqslant k \leqslant n} \big\vert \eta_k \big\vert \rightarrow 0~in~probability,\) and for every \(\epsilon > 0\), there exist \(n^*\) and \(\delta > 0\) such that for all \(n \geqslant n^*\), \({\mathsf{P}}\left\{ \max_{1 \leqslant k \leqslant n\delta} |\eta_{n+k} - \eta_n| > \epsilon \right\} < \epsilon.\)

Let \(A =\{A(t;\lambda), \lambda \in \Lambda\}\) be a family of boundary functions for some index set \(\Lambda\). For each \(\lambda \in \Lambda\), define \[\begin{align} T = T_{\lambda} = \inf\{n \geqslant 1: Z_n > A(n;\lambda) \},~~~\inf \varnothing = \infty. \label{1462b} \end{align}\tag{41}\] It is easy to see that under the positive drift assumption \(\mu > 0\), we have \(T_\lambda < \infty\) for all \(\lambda > 0\) with probability one. In nonlinear renewal theory, one focuses on asymptotic approximations for the distribution of the overshoot and for the expected stopping time \({\mathsf{E}}[T]\) as the boundary tends to infinity.

To illustrate the motivation of investigating stopping times \(T_\lambda\) in (41 ), we consider the following simple example: let \(\eta_n = 0\) and \(A(n;\lambda)= \lambda n^\alpha\) for \(0 \leqslant\alpha < 1\), then \[\begin{align} \label{1463a} T:=T_\lambda = \inf\{n \geqslant 1: S_n > \lambda n^\alpha \}. \end{align}\tag{42}\] Note that (42 ) is a standard formulation in nonlinear renewal theory, developed in [47], [48], which was originally motivated by problems arising in sequential analysis for statistical models. This formulation plays a central role in analyzing boundary crossing probabilities and provides the foundation for deriving asymptotic approximations of stopping times in a wide range of statistical applications.

By combining renewal theory with Wald’s identity, a standard approach is to investigate the difference between \(T_\lambda\) and a stopping time defined by crossing linear boundaries with varying drift. Specifically, we define \[\label{1463b} \tau(c,u) = \inf\{n \geqslant 1 : S_n - nu > c\}, \quad c \geqslant 0,~ 0 < u \leqslant\mu,\tag{43}\] and aim to establish the uniform integrability of \(|T_{\lambda} - \tau(c_{\lambda},d_{\lambda})|^p\), \(p \geqslant 1\), for suitable choices of \(c_{\lambda}\) and \(d_{\lambda}\). Nonlinear renewal theory is then derived directly from the corresponding results in the linear case with varying drift, by leveraging uniform integrability and the weak convergence of the overshoot.

When \(n\) is the first time of \(S_n\) crossing the boundary \(\lambda n^\alpha\), by Wald’s identity, we have \(n \mu = S_n \approx \lambda n^\alpha \Rightarrow n \approx \big(\lambda/\mu \big)^{1/(1-\alpha)}.\) Denote \(b_\lambda= \big(\lambda/\mu\big)^{1/(1-\alpha)}\). By using linearlization, \({d( \lambda n^\alpha)}/{d \alpha} = \alpha \lambda n^{\alpha - 1} \approx \alpha \mu,\) we can approximate the curve boundary \(\lambda n^\alpha\) by \(l_\lambda(n)= \mu b_\lambda + \alpha \mu (n - b_\lambda),\) Taylor’s expansion at \(b_\lambda\). Therefore, \(c = c_\lambda = \mu (1- \alpha) b_\lambda > 0\) and \(d = d_\lambda = \alpha \mu\) in this simpe case. Since we need to find upper and lower bounds of the curve boundary for each \(n\), \(\alpha\) is in a range of \([0, 1)\). This implies that we consider a random walk with varying drift or a uniform renewal theory. The reader is referred to the above mentioned articles for details.

Next, we present a multivariate nonlinear renewal theory from [54] as follows: let \(\{(X_n, Y_n),n=1,2,\cdots\}\) be a sequence of i.i.d.random variables in \({\boldsymbol{R}}^{d+1}\), where \(X_1 \in {\boldsymbol{R}^1}\) with positive drift and finite variance as before, and \(Y_1 \in {\boldsymbol{R}}^d\) with \({\mathsf{E}}[Y_1]=\vec{0}\) and \(\mathit{Var}(Y_1) =\Sigma_Y\). Denote \(\Sigma\) as the variance–covariance matrix of \((X_1, Y_1)\). Let \(S_n = \sum_{k=1}^n X_k\), and \({\boldsymbol{W}}_n = \sum_{k=1}^n Y_k\) with \({\boldsymbol{W}}_0=\vec{0}\).

Consider the stopping time \[\label{tau95general} \tau = \tau_b := \inf \{ n \geqslant 0: S_n - H({\boldsymbol{W}}_n+n\epsilon_n) > b\},\tag{44}\] where \(\epsilon_n \in {\boldsymbol{R}}^d\) with \(\epsilon_n \rightarrow 0_d\) as \(n \to \infty\), and \(H: {\boldsymbol{R}}^d \rightarrow {\boldsymbol{R}}\) with \(H(\vec{0}) = 0\). The \(\epsilon_n\) captures potential additional perturbations.

To simplify the presentation, we refer the reader to [54] for the complete version of these renewal theory results, and restrict ourselves here to the necessary special cases. The first result states that, under suitable regularity conditions and normalization, \(({\boldsymbol{W}}_\tau,\tau)\) is asymptotically distributed as a \((d+1)\)-dimensional normal random vector. Intuitively, this asymptotic normality arises from a renewal-type central limit theorem, where the cumulative effect of many small increments leads to a Gaussian limit. To describe this result more precisely, we introduce the following notation for characterizing the covariance matrix of \(({\boldsymbol{W}}_\tau,\tau)\). For any \(\nu \in \mathbf{R}^d\), define \[\begin{align} \label{notation-thm-prob} \notag M(\nu) = \left( \begin{matrix} 0 & I_d \\ -\mu_X^{-1} & \mu_X^{-1}\nu^t \end{matrix} \right), & ~~~ \tilde{\Sigma}^*(\nu) = M(\nu) \Sigma M(\nu)^t. \end{align}\tag{45}\] Here \(~^t\) denotes transpose. Further let \(\mathcal{N}_{b,\nu}\) follow a \((d+1)\)-dimensional normal distribution with mean \((0, \cdots, 0, {\mathsf{E}}[\tau_b])^t\) and covariance matrix \(\tilde{\Sigma}^*(\nu){\mathsf{E}}[\tau_b]\).

The characterization for \(\mathcal{N}_{b,\nu}\) in Theorem 7 requires the characterization of \({\mathsf{E}}[\tau_b]\), which is covered by the following theorem.

The results remain valid when \((X_1, Y_1)\) satisfies conditions (C1) and (C2), and when \(\epsilon_n\) satisfies condition (C5) in the Appendix of [54]; see Theorem A.1 therein for details. These results also extend to more general \(H\)-functions and yield the asymptotic distribution of the overshoot \(S_\tau - H({\boldsymbol{W}}_\tau + \tau \epsilon_\tau) - b.\) In particular, once the overshoot distribution is available, the constants \(\rho_+\) and \(c_0\) in Theorem 8 become explicitly computable (see the Appendix of [54]).

Last, we summarize a nonlinear Markov renewal theory from [55]. For convenience of notation, in the remaining part of this and the next subsection, we let \(\{(X_n, S_n),~ n \geqslant 0\}\) denote a Markov random walk on \({\cal X} \times \mathbf{R}\). That is, let \(\{X_n, n\geqslant 0\}\) be a Markov chain on a general state space \({\cal X}\) with \(\sigma\)-algebra \(\cal A\), which is irreducible with respect to a maximal irreducibility measure on \(({\cal X},\cal A)\) and is aperiodic. Let \(S_n = \sum_{k=1}^n \xi_k\) be the additive component, taking values on the real line \({\boldsymbol{R}}\), such that \(\{(X_n,S_n), n\geqslant 0\}\) is a Markov chain on \({\cal X} \times {\boldsymbol{R}}\) with transition probability \[\begin{align} \label{3461} & & {\mathsf{P}}\{(X_{n+1},S_{n+1}) \in A \times (B+s) | (X_n,S_n) = (x,s)\} \\ &=& {\mathsf{P}}\{(X_1,S_1) \in A \times B | (X_0,S_0) = (x,0)\} = {\mathsf{P}}(x,A \times B), \nonumber \end{align}\tag{48}\] for all \(x \in {\cal X},~ A \in {\cal A}\) and \(B \in {\cal B}({\boldsymbol{R}})\) (:= Borel \(\sigma\)-algebra on \({\boldsymbol{R}}\)). The chain \(\{(X_n,S_n), n \geqslant 0 \}\) is called a Markov random walk. In this subsection, let \({\mathsf{P}}_\nu~({\mathsf{E}}_{\nu})\) denote the probability (expectation) under the initial distribution on \(X_0\) being \(\nu\). If \(\nu\) is degenerate at \(x\), we shall simply write \({\mathsf{P}}_x~({\mathsf{E}}_x)\) instead of \({\mathsf{P}}_\nu ~({\mathsf{E}}_{\nu})\). We assume throughout this subsection that there exists a stationary probability distribution \(\pi\), \(\pi(A) = \int {\mathsf{P}}(x,A) {\mathrm{d}}\pi(x)\) for all \(A \in {\cal A}\) and \({\mathsf{E}}_{\pi} [\xi_1] >0\).

Let \(\{Z_n = S_n + \eta_n, n\geqslant 0\}\) be a perturbed Markov random walk in the following sense: \(S_n\) is a Markov random walk, \(\eta_n\) is \({\cal F}_n\)-measurable, where \({\cal F}_n\) is the \(\sigma\)-algebra generated by \(\{(X_k,S_k), 0\leqslant k \leqslant n\}\), and \(\eta_n\) is slowly changing. Let \(\{A =A(t;\lambda), \lambda \in \Lambda\}\) be a family of boundary functions for some index set \(\Lambda\). Define \[\begin{align} T = T_{\lambda} = \inf\{n \geqslant 1: Z_n > A(n;\lambda) \},~~~\inf \varnothing = \infty,~for~each~\lambda \in \Lambda. \end{align}\] It is easy to see that for all \(\lambda > 0\), \(T_\lambda < \infty\) with probability \(1\). This section concerns the approximations of the distribution of the overshoot and expected stopping time \({\mathsf{E}}_\nu [T]\) as the boundary tends to infinity.

A Markov chain \(\{X_n,n \geqslant 0\}\) on a state space \({\cal X}\) is called \(V\)-uniformly ergodic if there exists a measurable function \(V: {\cal X} \rightarrow [1,\infty)\), with \(\int V(x){\mathrm{d}}\pi(x) < \infty\), such that, for any Borel measurable function \(h\) on \({\cal X}\) satisfying \(||h||_V:= \sup_x |h(x)|/V(x) < \infty\), we have \[\begin{align} &~& \lim_{n \rightarrow \infty} \sup_{x \in {\cal X}} \bigg\{ \frac{|{\mathsf{E}}[h(X_n)|X_0=x] - \int h(x){\mathrm{d}}\pi(x)|}{V(x)} : x \in {\cal X}, |h| \leqslant V \bigg\} = 0. \label{ue} \end{align}\tag{49}\] In this subsection, we shall assume that \(\{X_n,n \geqslant 0\}\) is \(V\)-uniformly ergodic. Under irreducibility and aperiodicity assumption, \(V\)-uniform ergodicity implies that there exist \(r> 0\) and \(0 < \rho <1\) such that for all \(h\) and \(n \geqslant 1\), \[\begin{align} \sup_{x \in {\cal X}} \frac{|{\mathsf{E}}[h(X_n)|X_0=x ]- \int h(y) {\mathrm{d}}\pi(y)|}{V(x)} \leqslant r {\rho}^n \|h\|_V; \end{align}\] see pages 382-383 of [56]. When \(V \equiv 1\), this reduces to the classical uniform ergodicity condition.

The following assumptions for Markov chains will be used in this subsection.

A1. \(\sup_x \big\{\frac{{\mathsf{E}}[V(X_1)]}{V(x)} \big\} < \infty\),

A2. \(\sup_x {\mathsf{E}}_x [|\xi_1|^2] < \infty\) and \(\sup_x \big\{\frac{{\mathsf{E}}[|\xi_1|^r V(X_1)]}{V(x)} \big\} < \infty\) for some \(r \geqslant 1\).

A3. Let \(\nu\) be an initial distribution of the Markov chain \(\{X_n, n \geqslant 0\}\), assume that for some \(r \geqslant 1\), \[\sup_{||h||_V \leqslant 1} |\int_{x \in {\cal X}} h(x) {\mathsf{E}}_x[|\xi_1|^{r}] {\mathrm{d}}\nu(x) | < \infty.\]

A Markov random walk is called lattice with span \(d > 0\) if \(d\) is the maximal number for which there exists a measurable function \(\gamma: {\cal X} \to [0,\infty)\) called the shift function, such that \({\mathsf{P}}\{\xi_1-\gamma(x) + \gamma(y) \in \{\cdots, -2d, -d, 0, d, 2d,\cdots\}| X_0 = x, X_1 = y\} = 1\) for almost all \(x,y\in {\cal X}\). If no such \(d\) exists, the Markov random walk is called nonlattice. A lattice random walk whose shift function \(\gamma\) is identically 0 is called arithmetic.

To establish the nonlinear Markov renewal theorem, we shall make use of (48 ) in conjunction with the following extension of Cramer’s (strongly nonlattice) condition: there exists \(\delta > 0\) such that for all \(m,n=1,2,\cdots\), \(\delta^{-1} <m< n\), and all \(\theta \in R\) with \(|\theta| \geqslant\delta\) \[\begin{align} {\mathsf{E}}_{\pi} | {\mathsf{E}}[\exp(i\theta (\xi_{n-m}+\cdots+\xi_{n+m}))|X_{n-m},\cdots,X_{n-1},X_{n+1},\cdots,X_{n+m},X_{n+m+1}]| \leqslant e^{-\delta}. \end{align}\]

By using Markov renewal theory [57], [58] together with Wald’s equations for Markov random walks [59], our approach investigates the difference between \(T_{\lambda}\) and a stopping time defined by crossing linear boundaries with varying drift. Specifically, we define \[\label{tau-def} \tau(c,u) = \inf\{n \geqslant 1 : S_n - un > c\}, \quad c \geqslant 0,~ u \leqslant{\mathsf{E}}_\pi [\xi_1],\tag{50}\] and establish the uniform integrability of \(|T_{\lambda} - \tau(c_{\lambda},d_{\lambda})|^p\) for \(p \geqslant 1\) for suitable choices of \(c_{\lambda}\) and \(d_{\lambda}\). Nonlinear Markov renewal theory is then derived directly from the corresponding results in the linear case, by leveraging uniform integrability and the weak convergence of the overshoot.

Let \({\mathsf{P}}_+^u(x,B \times R) = {\mathsf{P}}_x\{X_{\tau(0,u)} \in B \}\) for \(u \leqslant{\mathsf{E}}_\pi [\xi_1]\), denote the transition probability associated with the Markov random walk generated by the ascending ladder variable \(S_{\tau(0,u)}\). Under the \(V\)-uniform ergodicity condition and \({\mathsf{E}}_\pi [\xi_1]> 0\), a similar argument as on page 255 of [57] yields that the transition probability \({\mathsf{P}}_+^u(x, \cdot \times R)\) has an invariant measure \(\pi_+^u\). Let \({\mathsf{E}}_{+}^u\) denote expectation under \(X_0\) having the initial distribution \(\pi_+^u\). When \(u={\mathsf{E}}_\pi [\xi_1]\), we denote \({\mathsf{P}}_+^{E_\pi \xi_1}\) as \({\mathsf{P}}_+\), and \(\tau_+ = \tau(0,{\mathsf{E}}_\pi [\xi_1])\). Define \[\begin{align} b &=& b_\lambda = \sup\{t\geqslant 1: A(t,\lambda) \geqslant t {\mathsf{E}}_\pi [\xi_1]\},~~~ \sup\varnothing = 1, \\ d &=& d_\lambda = (\frac{\partial A}{\partial t} )(b_\lambda;\lambda), \\ \overline{d} &=& \sup\{(\frac{\partial A}{\partial t} )(t;\lambda);~ t\geqslant b_\lambda,~ \lambda\in\Lambda\}, \\ R &=& R_\lambda = Z_T - A(T; \lambda), \\ R(c,u) &=& S_{\tau(c,u)} - u\tau(c,u) - c,~~~u\leqslant{\mathsf{E}}_\pi [\xi_1],~ c\geqslant 0, \\ r(u) &=& {\mathsf{E}}_{+}^u [R^2(0,u)]/2{\mathsf{E}}_{_+}^{u} [R(0,u)],~~~ u\leqslant{\mathsf{E}}_\pi [\xi_1], \\ G(r,u) &=& \int_r^\infty {\mathsf{P}}_{+}^u \{R(0,u) > s\}{\mathrm{d}}s/ {\mathsf{E}}_{+}^u [R(0,u)],~~~ u\leqslant{\mathsf{E}}_\pi [\xi_1],~ r\geqslant 0. \end{align}\]

We shall assume that \(A(t;\lambda)\) is twice differentiable in \(t\) and \(b_\lambda\) is finite so that \(d\) and \(\overline{d}\) are well defined. The next theorem is a Blackwell-type nonlinear Markov renewal theorem.

To study uniform integrabilitiy of the powers of the differences for linear and nonlinear stopping times, we shall first give the regularity conditions on \(\eta= \{\eta_n, n\geqslant 1\}\). The process \(\eta\) is said to be regular with \(p\geqslant 0\) and \(1/2 < \alpha \leqslant 1\) if there exists a random variable \(L\), a function \(f(\cdot)\) and a sequence of random variables \(U_n,~ n\geqslant 1\), such that \[\begin{align} &~& \eta_n = f(n)+U_n, ~ for~ n\geqslant L~and~ \sup_{x \in {\cal X}} {\mathsf{E}}_x [L^p] <\infty, \\ &~& \max_{1\leqslant j\leqslant\sqrt n} |f(n+j) - f(n)| \leqslant K,~~~ K < \infty, \\ &~& \big\{ \max_{1\leqslant j\leqslant n^\alpha} |U_{n+j}|^p,~ n\geqslant 1 \big\}~ is~uniformly~integrable, \\ &~& n^p \sup_{x \in {\cal X}} {\mathsf{P}}_x \big \{\max_{0\leqslant j\leqslant n} U_{n+j} \geqslant\theta n^\alpha \big\} \to 0~~ as~ n\to \infty,~{\rm for~all}~\theta > 0, \end{align}\] and for some \(w > 0,~ w < {\mathsf{E}}_\pi [\xi_1] - \overline{d}~ if~\alpha = 1\), \[\begin{align} \sum_{n=1}^\infty n^{p-1}\sup_{x \in {\cal X}} {\mathsf{P}}_x\{ -U_n \geqslant w n^\alpha\} < \infty. \end{align}\]

We shall set \(f(n)\) to be the median of \(\eta_n\) when \(\eta\) is not regular and extend \(f\) to a function on \([1,\infty)\) by linear interpolation. Therefore we can define \(\tau = \tau_\lambda = \tau(c_\lambda, d_\lambda)\) and \(c_\lambda = b_\lambda({\mathsf{E}}_\pi [\xi_1] - d_\lambda) - f(b_\lambda).\)

We need the following notations and definitions before Theorem 11.

For a given Markov random walk \(\{(X_n,S_n), n \geqslant 0\}\), let \(\nu\) be an initial distribution of \(X_0\), and define \[\nu^*(B) = \sum_{n=0}^\infty {\mathsf{P}}_\nu(X_n \in B), \quad B \in \mathcal{A}.\]

Let \(g = {\mathsf{E}}[\xi_1 \mid X_0, X_1]\) with \({\mathsf{E}}_\pi[|g|] < \infty\). Define operators \({\boldsymbol{P}}\) and \({\boldsymbol{P}}_\pi\) by \[({\boldsymbol{P}} g)(x) = {\mathsf{E}}_x[g(x,X_1,\xi_1)], \qquad {\boldsymbol{P}}_\pi g = {\mathsf{E}}_\pi[g(X_0,X_1,\xi_1)],\] and set \(\overline{g} = {\boldsymbol{P}} g\).

We shall consider solutions \(\Delta(x) = \Delta(x; g)\) of the Poisson equation \[\label{34628} \big(I - {\boldsymbol{P}}\big) \Delta = \big(I - {\boldsymbol{P}}_\pi \big) \overline{g}, \quad \nu^*\text{-a.s.}, \qquad {\boldsymbol{P}}_\pi \Delta = 0,\tag{57}\] where \(I\) is the identity operator.

Under conditions A1–A4, it is known [56] that the solution \(\Delta\) of 57 exists and is bounded.

When \(A(t,\lambda)= \lambda\), we have

5.2 Multivariate Markov Renewal Theory↩︎

In this subsection, we consider the multivariate Markov renewal theory developed in [57]. For a comprehensive background on multivariate renewal theory, the reader is referred to that paper and the references therein. We begin by presenting multivariate Markov renewal theorems with convergence rates. Recall that \(\{X_n, n \geqslant 0\}\) is assumed to be irreducible (with respect to some measure on \({\cal A}\)), aperiodic, and \(V\)-uniformly ergodic, and \(\{(X_n, S_n := \sum_{k=1}^n \xi_k),~ n \geqslant 0\}\) is a Markov random walk on \({\cal X} \times \mathbf{R}^d\). In addition, we assume \[\begin{align} &~& \sup_x\{\frac{{\mathsf{E}}_x[V(X_1)]}{V(x)} \} < \infty, \tag{58} \\ &~& \sup_x {\mathsf{E}}_x[|\xi_1|^2] < \infty~{\rm and~} \sup_x \Big\{ \frac{{\mathsf{E}}_x[|\xi_1|^r V(X_1)]}{V(x)} \Big\} < \infty. ~~{\rm for~ some}~ r \geqslant 2. \tag{59} \end{align}\]

Let \(\pi\) be the stationary distribution of \(\{X_n, n \geqslant 0\}\) and let \({\mathsf{P}}_\pi\) denote \(\int {\mathsf{P}}_x {\mathrm{d}}\pi(x)\) and \({\mathsf{E}}_\pi\) be expectation under \({\mathsf{P}}_\pi\). Hereafter, we use column vectors to denote \(\theta \in {\boldsymbol{R}}^d\), \(\theta^t\) to denote the transpose of \(\theta\), and \(|\theta|\) to denote its Euclidean norm \((\theta^t\theta)^{1/2}\).

Let \(\mu = {\mathsf{E}}_\pi[ \xi_1]\) and \(\Sigma = \lim_{n \to \infty} n^{-1} {\mathsf{E}}_\pi [\{(S_n-n\mu)(S_n- n\mu)^t\}]\), which are well defined under (49 ), 58 and (59 ). Let \(S_{n,j}\) (or \(\xi_{n,j}\), \(\mu_j\), \(\theta_j\)) denote the \(j\)th component of the \(d\)-dimensional vector \(S_{n}\) (or \(\xi_{n}\), \(\mu\), \(\theta\)). Suppose \(\mu_1 > 0\). Without loss of generality, it will be assumed that \(\Sigma\) is positive definite (i.e., \(\xi_n\) is strictly \(d\)-dimensional under \(\pi\)), because otherwise we can consider a lower-dimensional subspace instead. In the case \(d > 1\) define \[\gamma ={\mathsf{E}}_\pi[\{(\xi_{n,2}/\mu_1,\cdots,\xi_{n,d}/\mu_1)^t\}],~~~ \tilde{\Sigma} = (-\gamma, \;I_{d-1})\Sigma \left(\begin{array}{c} -\gamma^t \\I_{d-1} \end{array} \right),\] where \(I_k\) is the \(k\times k\) identity matrix. Note that \(\tilde{V}\) is the asymptotic covariance matrix (under \({\mathsf{P}}_\pi\)) of \(\{(S_{n,2},\cdots,S_{n,d})^t - S_{n,1} \gamma\} / \sqrt n\). For \(s \in {\boldsymbol{R}}^d\), define \(\tilde{s}=(s_2,\cdots,s_d)^t - s_1\gamma\).

First consider the case of i.i.d. \(\xi_n\), with \(d > 1\) and \(S_0=0\). The renewal measure is defined by \(U(B)=\sum_{n=0}^{\infty} {\mathsf{P}}\{S_n\in B\},\) and multivariate renewal theory is concerned with approximating \(U(s + \cdot)\) by \(\Psi_k(s + \cdot)\) as \(s_1 \to \infty\), where \(\Psi_k\) is a \(\sigma\)-finite measure on \({\boldsymbol{R}}^d\) whose density function (i.e., Radon-Nikodym derivative) with respect to Lebesgue measure is of the form \[\begin{align} \label{psi} \psi_k(s) = \frac{1}{\mu_1 \sqrt{det \tilde{\Sigma}}} (\frac{\mu_1}{2 \pi s_1})^{(d-1)/2} e^{-\mu_1 \tilde{s}^t \tilde{\Sigma}^{-1} \tilde{s}/2s_1} \{1 + \sum_{j=1}^k s_1^{-j/2} \omega_j(\tilde{s}/\sqrt{s_1})\} \end{align}\tag{60}\] for \(s_1 > 0\), and \(\psi_k(s)=0\) for \(s_1 \leqslant 0\), where \(\omega_j(u) = \sum_{l=0}^{n_j} q_l(u)\) and \(q_l(u)\) is a polynomial of degree \(l\) in \(u\) whose coefficients are associated with the Taylor expansion of \((1- {\mathsf{E}}[ e^{i \theta^t \xi_1}])^{-1}\) near \(\theta = 0\). For Markov random walks, the renewal measure involves not only \(\{S_n\}\) but also \(\{X_n\}\). For \(A \in {\cal A}\) and \(B \in {\cal B}\), define \[\begin{align} \label{RM} U_\nu^A(B)= \sum_{n=0}^{\infty} {\mathsf{P}}_\nu\{X_n\in A,S_n\in B \}. \end{align}\tag{61}\] We can approximate \(U_\nu^A(s + \cdot)\) by \(\pi(A) \Psi_k^{A,\nu}(s + \cdot)\), in which \(\Psi_k^{A,\nu}\) is a \(\sigma\)-finite measure on \({\boldsymbol{R}}^d\) with density function \(\psi_k^{A,\nu}\) with respect to Lebesgue measure, where \(\psi_k^{A,\nu}(s)=0\) for \(s_1 \leqslant 0\) and \(\psi_k^{A,\nu}(s)\) is given by (60 ) for \(s_1 > 0\), with the coefficients of the polynomials \(\omega_1(\tilde{s}),\cdots, \omega_k(\tilde{s})\) depending also on \(A\) and \(\nu\) via Taylor’s expansion of the Fourier transform of \(U_\nu^A\) near the origin, assuming that \[\begin{align} \label{mc} {\mathsf{E}}_\nu[V(X_1)(1+|S_1|^r)] < \infty \end{align}\tag{62}\] for some sufficiently large \(r\) (depending on \(k\)). Note that when \(\nu\) is degenerate at \((x,0)\), (62 ) follows from (59 ). The precise definition of \(\omega_j\) is given in Section 4.1 in [57], where we also prove the following multidimensional Markov renewal theorem with bounds on the remainders in approximating \(U_\nu^A(s + \cdot)\) by \(\pi(A) \Psi_k^{A,\nu}(s + \cdot)\) as \(s_1 \to \infty\), recalling the assumption \(\mu_1 > 0\).

For \(\varepsilon > 0\) and \(f:{\boldsymbol{R}}^d \to {\boldsymbol{R}}\), define the oscillation function \(\Omega_f(s;\varepsilon) = \sup\{ |f(s)-f(t)|: |s-t| \leqslant\varepsilon\}.\) Let \({\cal F}_b\) be the set of all Borel functions \(f:{\boldsymbol{R}}^d \rightarrow [0,1]\) such that \(f(s) = 0\) whenever \(s_1 \not\in [b,b+h],\) with fixed \(h >0\).

It is known that renewal theorems are often applied to the ladder random walk. The techniques used by [56] to prove the \(V\)-uniform ergodicity of a rich class of time series and queuing models can also be applied to show that their ladder random walks indeed satisfy conditions (49 ), 58 , (59 ) and (62 ). Recall that the (positive) ladder epoch of a Markov random walk \(\{(X_n,S_n), n \geqslant 0\}\) taking values in \({\cal X} \times {\boldsymbol{R}}\) is defined by \[\begin{align} \label{ladder} \tau_+ = \inf\{n \geqslant 1:S_n > 0\}. \end{align}\tag{63}\] For \(A \in {\cal A}\) and Borel subset \(B\) of \((0,\infty)\), define \[\begin{align} \label{lmc} {\mathsf{P}}_+(x, A \times B) = {\mathsf{P}}\{X_{\tau_+} \in A, S_{\tau_+} \in B|X_0=x \}. \end{align}\tag{64}\] The kernel \({\mathsf{P}}_+\) is the transition probability kernel of a Markov random walk that has the ladder chain as the underlying Markov chain on \({\cal X}\).

Next, we consider limit theorems for first passage times of Markov random walks by making use of Markov renewal theory for the ladder random walk, with kernel (64 ) in which \(\tau_+\) is defined by (63 ) in the case \(d=1\) and by \(\tau_+ = \inf\{n \geqslant 1: S_{n,1} > 0\}\) in the case \(d > 1\). It is assumed throughout this section that \({\mathsf{P}}_x(\tau_+ < \infty) =1\) for all \(x \in {\cal X}\) and that the ladder random walk is strongly nonlattice and satisfies conditions (49 ), 58 , (59 ) and 62 . Let \(\pi_+\) denote the invariant measure of the kernel \({\mathsf{P}}_+(x, A \times {\boldsymbol{R}}^d)\) which is assumed to be irreducible and aperiodic. Let \(\tau_1 = \tau_+, \;\tau_{j+1} = \inf\{n > \tau_j : S_{n,1} > S_{\tau_{j,1}}\}\) and \[\begin{align} T_b &=& \inf\{n \geqslant 1:S_{n,1} > b \}, \label{TB} \\ \mu^* &=& {\mathsf{E}}_{\pi_+} [S_{\tau_+}],~\Sigma_+ = \lim_{n \rightarrow \infty} n^{-1} {\mathsf{E}}_{\pi_+} [(S_{\tau_n}- n \mu^*) (S_{\tau_n}- n \mu^*)^t]. \end{align}\tag{65}\] Define \(\gamma_+\) and \(\widehat{\Sigma}_+\) as in (2.6) but with \(\pi_+,S_{\tau_n} - S_{\tau_{n-1}}, \Sigma_+\) in place of \(\pi,\xi_n\) and \(\Sigma\).

We begin by analyzing the asymptotic distribution of \((X_{T_b}, S_{T_b})\) for a Markov random walk. Let \(\sigma_+^2 = \lim_{n \to \infty} n^{-1} {\mathsf{E}}_\pi [S_{\tau_n,1} - n\mu_1^*]^2\).

For applications to nonlinear first passage problems, we replace \(S_n\) in 65 and in Theorem 14 by \(R_n = S_n + \Delta_n\), where \(\Delta_n\) represents some nonlinear perturbation. For the case of i.i.d. increments \(\xi_n\) and \(d=1\), such extension has been developed by [47]. The following theorem extends their result to the Markov case and \(d \geqslant 1\).

We next describe an important class of examples that motivate Theorem 15. Let \((X_n,S^*_n)\) be a Markov random walk such that \(X_n\) has stationary distribution \(\pi\) and \(\mu_{\pi}={\mathsf{E}}_\pi [S_1^*] \in {\boldsymbol{R}}^k\). Suppose \(g:{\boldsymbol{R}}^k \to {\boldsymbol{R}}\) and \(h:{\boldsymbol{R}}^k \to {\boldsymbol{R}}^{d-1}\) are twice continuously differentiable in some neighborhood of \(\mu_\pi\). Let \(R_{n,1} = n g(S_n^*/n),~(R_{n,2},\cdots,R_{n,d})^t = n h(S_n^\ast/n).\) Then Taylor expansions of \(g\) and \(h\) show that \(R_n\) can be expressed as \(S_n + \Delta_n\), where \(S_{n,1} = Dg(\mu_\pi)(S_n^* - n\mu_\pi)\) and \((S_{n,2},\cdots,S_{n,d})^t = Dh(\mu_\pi)(S_n^* - n\mu_\pi)\), in which \(Dg=(\partial g/\partial s_1,\cdots,\partial g/\partial s_k)\) and \(Dh = (\partial h_i/\partial s_j)_{1 \leqslant i \leqslant d-1,1 \leqslant j \leqslant k}.\) See [47], [48] for certain special cases of \(g\) when \(S^*_n\) has i.i.d. increments.

Building on these approximations, [57] derive explicit asymptotic expansions for the distribution of \((X_{T_b}, S_{T_b})\) as the boundary \(b \to \infty\). The key idea is to treat the Markov random walk as a “locally linear” process near the boundary crossing, using the ladder process to capture the first-passage events. This approach allows the otherwise complex dependency structure to be expressed in terms of the invariant measure \(\pi_+\) and the polynomial corrections \(\omega_1(\tilde{s}/\sqrt{s_1}; x, x_0)\).

For the ladder random walk, the measure \(U_+\) and its density \(u_+(x;x_0,B)\) play a central role. Intuitively, \(u_+(x;x_0,B)\) describes the distribution of the state when the random walk first exceeds the boundary in the positive direction. The approximation \(u_+(x;x_0,B) \approx p_+(x) \, \psi_1^{x,x_0}(s)\) provides a tractable way to incorporate the effects of both the initial state \(x_0\) and the specific “overshoot” \(s\) over the boundary.

These expansions also make it possible to compute functionals of \((X_{T_b}, S_{T_b})\), such as moments or tail probabilities, with controlled error terms. In particular, the convergence rates of the polynomial corrections in \(\psi_1^{x,x_0}(s)\) give explicit bounds on the accuracy of the approximations, highlighting the role of the Markov dependence structure and the moment conditions on \(\xi_n\).

By expressing the first passage distribution in terms of the ladder process and its invariant measure, one can separate the long-term stationary behavior of the Markov chain from the short-term fluctuations of the random walk. This separation simplifies the analysis and allows classical renewal-theoretic techniques, such as expansions in powers of \(1/\sqrt{b}\), to be applied in the multivariate and Markov-dependent setting.

In summary, the multivariate Markov renewal framework provides a powerful tool for approximating the distribution of boundary crossing times and associated states, even in complex dependent settings. The combination of absolute continuity conditions, density approximations, and polynomial corrections offers a flexible and explicit method to handle both finite and general state spaces.

If one ignores terms of the order \(O(1/\sqrt{b})\) in the integral in Theorem 7, then the integral reduces via integration by parts to \[\begin{align} &~& {\mathsf{P}}\{W \in C\} \int_{\cal X} \int_{y \in A} \int_{s \in {\boldsymbol{R}}^d,s_1 > u} (s_1 -u) {\mathsf{P}}_+(x,{\mathrm{d}}y \times {\mathrm{d}}s){\mathrm{d}}\pi_+(x) \\ &=& {\mathsf{P}}\{W \in C\} \int_{u}^{\infty} {\mathsf{P}}_{\pi_+}\{X_{\tau_+} \in A,~ S_{\tau_+,1} > v\}{\mathrm{d}}v, \end{align}\] where \(W\) is a \((d-1)\)-dimensional Gaussian vector with mean \(0\) and covariance matrix \(\widehat{\Sigma}_+\), as in Theorem 14.

In addition to the multivariate Markov renewal theory developed in [57], [59] established Wald’s equations, as well as results on first-passage times and moments of ladder variables in Markov random walks. Uniform Markov renewal theory can be found in [58], [60]. Applications to sequential analysis are discussed in [9], while applications to changepoint detection in hidden Markov models appear in [55], [61], [62]. Extending these results to multivariate nonlinear Markov renewal theory remains an interesting and open problem.

6 SELECTED CONTRIBUTIONS TO BIOSTATISTICS↩︎

A major theme in the latter portion of Lai’s career was applying many of the tools described above from probability, sequential hypothesis testing, and changepoint detection to designing biomedical clinical trials and other biostatistics problems. In this section we review just a small subset of Lai’s biostatistics work, on optimal group sequential designs (Section 6.1) and on adaptive designs for mid-trial sample size re-estimation (Section 6.2). A more complete view of Lai’s work in these and related areas can be found in his book Sequential Experimentation in Clinical Trials: Design and Analysis [63], whose aim was to provide a bridge from statistical theory to biostatistics practice, as evidenced through the popular R software package sp23design [64] implementing the phase II/III and other designs covered in the book. Lai’s work in this area also had a sizeable institutional impact, and Lai co-led Stanford Medical School’s Biostatistics Core and founded the Center for Innovative Study Design (CISD), producing a stream of adaptive/personalized-medicine trial methodology and applications. Some of Lai’s major areas of contribution not covered here are response-adaptive randomization [65], biomarker-guided and adapative-enrichment designs [66], comparative-effectiveness and point-of-care (POC) trials [67], and phase I dose-finding trials [68], [69].

6.1 Group Sequential Testing for Clinical Trials↩︎

The optimality theory discussed in Section 2 for fully-sequential testing would require modification for widespread use in clinical trials, where the dominant methodology is group sequential testing in which groups of patients are analyzed. [70] modified the preceding theory for group sequential tests in a one-parameter exponential family 10 of density functions, for which Hoeffding’s lower bound can be expressed as \[\label{a13} {\mathsf{E}}_\theta(T)\geqslant-\zeta^{-1}\log(\alpha+\beta)-\left(\zeta^{-2}\sigma/2\right)\left\{(\sigma/4)^2 -\zeta\log(\alpha+\beta)\right\}^{1/2}+\zeta^{-2}\sigma^2/8\tag{66}\] where \(\sigma^2=(\theta_1-\theta_0)^2 b^{''}(\theta)=\mathop{\mathrm{\mathsf{Var}}}_\theta\{(\theta_1-\theta_0)X_i\},\; \zeta=\max\{I(\theta,\theta_0),I(\theta,\theta_1)\}\), and \[I(\theta,\lambda)={\mathsf{E}}_\theta\left[\log\{f_\theta(X_i)/f_\lambda(X_i)\}\right] =(\theta-\lambda) b'(\theta)-\left( b(\theta)- b(\lambda)\right)\] is the Kullback–Leibler information number. The lower bound 66 does not take into consideration the fact that \(T\) can assume only several possible values in the case of group sequential designs. The first step of [70] is to take this into consideration by providing an asymptotic lower bound for \(T\) in the following theorem. Let \(n_0=0\).

The \(n_j\) in Theorem 17 can in fact be random variables independent of \(X_1,X_2,\dots\). In this case the preceding argument can still be applied after conditioning on \((n_1,\dots,n_k)\). The next step of [70] is to extend Theorem 2 about asymptotic optimality of Lorden’s 2-SPRT to the group sequential setting in the following. Let \(S_n=\sum_{i=1}^n X_i\).

Whereas the theory for the group sequential 2-SPRT in Theorem 18 requires specification of \(\theta\), the group sequential GLR in the next section replaces \(\theta\) at the \(i\)th interim analysis by the maximum likelihood estimate \(\hat{\theta}_{n_i}\), similar to the tests in Section 2.3 for the fully sequential case. The resulting group sequential GLR test of [70] attains the asymptotic lower bound 69 at every fixed \(\theta\) and that its power is comparable to the upper bound \(1-\beta\) at \(\theta_1\), under the assumption that the group sizes satisfy 67 with \(n_k \sim |\log\alpha|/I(\theta^*,\theta_0)\), as \(\alpha+\beta\to0\) such that \(\log\alpha\sim\log\beta\).

6.2 Efficient Adaptive Designs for Mid-Course Sample Size Adjustment↩︎

In this section, we review some of Lai’s and his coathors’ work on re-estimating the sample size for a group sequential test midway through a clinical trial. There was increasing interest in this topic because of the ethical and economic considerations in the design of clinical trials to test the efficacy of new treatments and lack of information on the new treatments being tested. One challenge in this area is the difficulty in specifying the alternative on which the power of a test is based, and pilot studies may be unavailable or difficult to interpret, thus there was interest in finding ways to adjust the test statistics while maintaining control on the the type I error probability.

Much of the existing literature focused on 2-stage methods for testing about the normal mean. A unified treatment, developed by [71], [72] in the general framework of multiparameter exponential families uses efficient GLR statistics and adds a third stage to adjust for the sampling variability of the first-stage parameter estimates that determine the second-stage sample size. The possibility of adding a third stage to improve two-stage designs dated back to [73]. Whereas Lorden’s upper bounds for the type I error probability are too conservative for clinical trial applications which must follow strict regulatory specifications, [71] overcame this difficulty by modifying the numerical methods to compute the type I error probability, and also extended the three-stage test to multiparameter and multi-armed settings, thus greatly broadening the scope of these efficient adaptive designs.

6.2.1 An Adaptive 3-stage GLR Test↩︎

[71], [72] consider the general framework of the multiparameter exponential family \(f_{\boldsymbol{\theta}}(x)=\exp(\boldsymbol{\theta}^T x-\psi(\boldsymbol{\theta}))\) considered in 24 . Let \(\Lambda_{i,j}\), \(j\in\{0,1\}\), denote the GLR statistic in this family comparing \(\boldsymbol{\widehat{\theta}}_{n_i}\) to the composite \(u(\boldsymbol{\theta})=u_j\) and \(u\) is a smooth real-valued function defining the hypotheses, below. Rather than considering simple null and alternative hypotheses, [71] use the GLR statistics \(\Lambda_{i,0}\) and \(\Lambda_{i,1}\) in an adaptive three-stage test of the composite null hypothesis \(H_0: u(\boldsymbol{\theta})\leqslant u_0\), such that \[\text{I(\boldsymbol{\theta},\boldsymbol{\lambda}) is increasing in |u(\boldsymbol{\lambda})-u(\boldsymbol{\theta})| for every fixed \boldsymbol{\theta}.}\] Let \(n_1=m\) be the sample size of the first stage and \(n_3=M\) be the maximum total sample size, both specified before the trial. Let \(u_1>u_0\) be the alternative implied by the maximum sample size \(M\) and the reference type II error probability \(\tilde{\alpha}\). That is, \(u_1(>u_0)\) is the alternative where the fixed sample size (FSS) GLR test with type I error probability \(\alpha\) and sample size \(M\) has power \(\inf_{\boldsymbol{\theta}:u(\boldsymbol{\theta})=u_1} {\mathsf{P}}_{\boldsymbol{\theta}}\{\text{Reject H_0}\}\) equal to \(1-\tilde{\alpha}\). The three-stage test of \(H_0: u(\boldsymbol{\theta})\leqslant u_0\) stops and rejects \(H_0\) at stage \(i\leqslant 2\) if \[\label{eq:6} n_i<M,\quad u\big(\boldsymbol{\widehat{\theta}}_{n_i}\big)>u_0,\quad\text{and}\quad\Lambda_{i,0}\geqslant b.\tag{70}\] Early stopping for futility (accepting \(H_0\)) can also occur at stage \(i\leqslant 2\) if \[\label{eq:7} n_i<M,\quad u\big(\boldsymbol{\widehat{\theta}}_{n_i}\big)<u_1,\quad\text{and}\quad\Lambda_{i, 1} \geqslant\tilde{b}.\tag{71}\] The test rejects \(H_0\) at stage \(i=2\) or \(3\) if \[\label{eq:8} n_i=M,\quad u\big(\boldsymbol{\widehat{\theta}}_M\big)>u_0,\quad\text{and}\quad\Lambda_{i,0}\geqslant c,\tag{72}\] accepting \(H_0\) otherwise. The sample size \(n_2\) of the three-stage test is given by \[n_2=m\vee\left\{M \wedge\left\lceil (1 + \rho_m)\, n\big(\boldsymbol{\widehat{\theta}}_m\big)\right\rceil\right\},\] with \[\label{eq:5} n(\boldsymbol{\theta}) = \min\left\{|\log\alpha|\Big/\inf_{\boldsymbol{\lambda}:u(\boldsymbol{\lambda})=u_0} I(\boldsymbol{\theta},\boldsymbol{\lambda}),\; |\log\tilde{\alpha}|\Big/\inf_{\boldsymbol{\lambda}:u(\boldsymbol{\lambda})=u_1}I(\boldsymbol{\theta},\boldsymbol{\lambda})\right\},\tag{73}\] where \(I(\boldsymbol{\theta},\boldsymbol{\lambda})\) is the Kullback–Leibler information number and \(\rho_m>0\) is an inflation factor to adjust for uncertainty in \(\boldsymbol{\widehat{\theta}}_m\). Note that 73 is an asymptotic approximation to Hoeffding’s lower bound 66 . Letting \(0<\varepsilon,\tilde{\varepsilon}<1\), define the thresholds \(b, \tilde{b}\), and \(c\) to satisfy the equations \[\begin{align} &\sup_{\boldsymbol{\theta}:u(\boldsymbol{\theta})=u_1}{\mathsf{P}}_{\boldsymbol{\theta}}\{\text{\eqref{eq:7} occurs for i=1 or 2}\}=\tilde{\varepsilon}\tilde{\alpha},\tag{74}\\ &\sup_{\boldsymbol{\theta}:u(\boldsymbol{\theta})=u_0}{\mathsf{P}}_{\boldsymbol{\theta}}\{\text{\eqref{eq:7} does not occur for i\leqslant 2, \eqref{eq:6} occurs for i=1 or 2}\}=\varepsilon\alpha,\\ &\sup_{\boldsymbol{\theta}:u(\boldsymbol{\theta})=u_0}{\mathsf{P}}_{\boldsymbol{\theta}}\{\text{\eqref{eq:6} and \eqref{eq:7} do not occur for i\leqslant 2, \eqref{eq:8} occurs}\}=(1-\varepsilon)\alpha.\tag{75} \end{align}\] The probabilities in 74 –75 can be computed by using the normal approximation to the signed-root likelihood ratio statistic \[\ell_{i,j}=\left\{\text{sign}\left(u(\boldsymbol{\widehat{\theta}}_{n_i})-u_j\right)\right\} (2n_i\Lambda_{i,j})^{1/2}\] (\(1\leqslant i\leqslant 3\) and \(j=0,1\)) under \(u(\boldsymbol{\theta})=u_j\). When \(u(\boldsymbol{\theta})=u_j\), \(\ell_{i,j}\) is approximately normal with mean 0, variance \(n_i\), and the increments \(\ell_{i,j}-\ell_{i-1,j}\) are asymptotically independent. We can therefore approximate \(\ell_{i,j}\) by a sum of independent standard normal random variables under \(u(\theta)=u_j\) and thereby determine \(b, \tilde{b}\), and \(c\). Note that this normal approximation can also be used for the choice of \(u_1\) implied by \(M\) and \(\tilde{\alpha}\).

A special multiparameter case of particular interest in clinical trials involves \(K\) independent populations having density functions \(\exp\{\theta_k x-\tilde{\psi}_k(\theta_k)\}\) so that \(\boldsymbol{\theta}^T \boldsymbol{x}-\psi(\boldsymbol{\theta})=\sum_{k=1}^K \{\theta_k x_k-\tilde{\psi}(\theta_k)\}\). In multi-arm trials, for which different numbers of patients are assigned to different treatments, the GLR statistic \(\Lambda_{i,j}\) for testing the hypothesis \(u(\theta_1,\dots,\theta_K)=u_j\) (\(j=0\) or \(1\)) at stage \(i\) has the form \[\Lambda_{i,j}=\sum_{k=1}^Kn_{ki}\left\{\widehat{\theta}_{k,n_{ki}}\bar{X}_{k,n_{ki}}- \tilde{\psi}\left(\widehat{\theta}_{k,n_{ki}}\right)\right\} -\sup_{\boldsymbol{\theta}: u(\theta_1,\dots,\theta_K)=u_j} \,\sum_{k=1}^K n_{ki} \left\{\theta_k\bar{X}_{k,n_{ki}}-\tilde{\psi}(\theta_k)\right\},\] in which \(n_{ki}\) is the total number of observations from the \(k\)th population up to stage \(i\). Letting \(n_i=\sum_{k=1}^K n_{ki}\), the normal approximation to the signed root likelihood ratio statistic is still applicable when \(n_{ki}=p_kn_i+O_p(n_i^{1/2})\), where \(p_1,\dots,p_K\) are nonnegative constants that sum up to 1, as in random allocation of patients to the \(K\) treatments (for which \(p_k=1/K\)); see [70].

6.2.2 Mid-Course Modification of Maximum Sample Size↩︎

The adaptive designs in the preceding section can be modified to accommodate the possibility of mid-course increase of the maximum sample size from \(M\) to \(\widetilde{M}\). Let \(u_2\) be the alternative implied by \(\widetilde{M}\) so that the level-\(\alpha\) GLR test with sample size \(\widetilde{M}\) has power \(1-\tilde{\alpha}\). Note that \(u_1>u_2>u_0\). Whereas the sample size \(n_3\) is chosen to be \(M\) in Section 6.2.1, we now define \[\begin{align} \tilde{n}(\boldsymbol{\theta})&=\min\left\{|\log\alpha|\Big/\inf_{\boldsymbol{\lambda}:u(\boldsymbol{\lambda})=u_0}I(\boldsymbol{\theta},\boldsymbol{\lambda}),\;|\log\tilde{\alpha}|\Big/\inf_{\boldsymbol{\lambda}:u(\boldsymbol{\lambda})=u_2}I(\boldsymbol{\theta},\boldsymbol{\lambda})\right\},\\ n_3&=n_2\vee\left\{M'\wedge\left\lceil(1+\rho_m)\tilde{n}\big(\hat{\boldsymbol{\theta}}_{n_2}\big)\right\rceil\right\}, \end{align}\] where \(M<M' \leqslant\widetilde{M}\) and \(n_2=m\vee \{M\wedge (1+\rho_m)\tilde{n}(\boldsymbol{\widehat{\theta}}_m)\}\). We can regard the test as a group sequential test with 4 groups and \(n_1=m\), \(n_4=\widetilde{M}\), but with adaptively chosen \(n_2\) and \(n_3\). If the test does not end at the third stage, continue to the fourth and final stage with sample size \(n_4=\widetilde{M}\). Its rejection and futility boundaries are similar to those in Section 6.2.1. Extending our notation \(\Lambda_{i,j}\) to \(1\leqslant i\leqslant 4\) and \(0\leqslant j\leqslant 2\), the test stops at stage \(i\leqslant 3\) and rejects \(H_0\) if \[n_i<\widetilde{M},\quad u\big(\boldsymbol{\widehat{\theta}}_{n_i}\big)>u_0,\quad\text{and}\quad\Lambda_{i,0}\geqslant b,\] stops and accepts \(H_0\) if \[\label{eq:15} n_i<\widetilde{M},\quad u\big(\hat{\boldsymbol{\theta}}_{n_i}\big)<u_2,\quad\text{and}\quad\Lambda_{i,2}\geqslant\tilde{b},\tag{76}\] and rejects \(H_0\) at stage \(i=3\) or \(4\) if \[n_i=\widetilde{M},\quad u\big(\hat{\boldsymbol{\theta}}_{\widetilde{M}}\big)>u_0,\quad\text{and}\quad \Lambda_{i,0}\geqslant c,\] accepting \(H_0\) otherwise. The thresholds \(b, \tilde{b}\), and \(c\) can be defined by equations similar to 74 –75 to insure the overall type I error probability to be \(\alpha\). For example, in place of 74 , \[\label{eq:16} \sup_{\boldsymbol{\theta}:u(\boldsymbol{\theta})=u_2} {\mathsf{P}}_\theta \{\text{\eqref{eq:15} occurs for some i\leqslant 3}\}=\tilde{\varepsilon}\tilde{\alpha}.\tag{77}\] The basic idea underlying 77 is to control the type II error probability at \(u_2\) so that the test does not lose much power there in comparison with the GLR test that has sample size \(\widetilde{M}\) (and therefore power \(1-\tilde{\alpha}\) at \(u_2\)).

6.2.3 Asymptotic Theory for the Adaptive 3-Stage GLR Test↩︎

[71], [72] established the asymptotic optimality of the above three-stage test in the following.

Since \(m\sim a|\log\alpha|\) and \(M\sim A|\log\alpha|\) and since the thresholds \(b\), \(\tilde{b}\), and \(c\) are defined by 74 –75 , [72] use an argument similar to the proof of Theorem 2(ii) of [70] to show that 79 holds. They then use the following argument to prove 80 . Let \(\Theta_0=\{\boldsymbol{\theta}:u(\boldsymbol{\theta})\leqslant u_0\}\), \(\Theta_1=\{\boldsymbol{\theta}:u(\boldsymbol{\theta})\geqslant u_1\}\). For \(i=0,1\), \[\inf_{\boldsymbol{\lambda} \in \Theta_i}I(\boldsymbol{\theta},\boldsymbol{\lambda})=I_i(\boldsymbol{\theta}), \qquad\text{where } I_i(\boldsymbol{\theta})=\inf_{\boldsymbol{\lambda}:u(\boldsymbol{\lambda})=u_i} I(\boldsymbol{\theta},\boldsymbol{\lambda}). \label{4464}\tag{81}\] Take any \(\boldsymbol{\lambda}\in\Theta_0\) and \(\tilde{\boldsymbol{\lambda}}\in\Theta_1\). From 81 and Hoeffding’s lower bound 66 , it follows that for a test that has error probabilities \(\alpha\) and \(\tilde{\alpha}\) at \(\boldsymbol{\lambda}\) and \(\tilde{\boldsymbol{\lambda}}\) and take at least \(m\) and at most \(M\) observations, its sample size \(T\) satisfies \[{\mathsf{E}}_{\boldsymbol{\theta}}(T)\geqslant m\vee\left\{M\wedge\frac{[1+o(1)]|\log\alpha|}{I_0(\boldsymbol{\theta})\vee I_1(\boldsymbol{\theta})}\right\} \label{4465}\tag{82}\] as \(\alpha+\tilde{\alpha}\to0\) such that \(\log\alpha\sim\log\tilde{\alpha}\). The second-stage sample size of the adaptive test is a slight inflation of the Hoeffding-type lower bound 82 with \(\boldsymbol{\theta}\) replaced by the maximum likelihood estimate \(\hat{\boldsymbol{\theta}}_m\) at the end of the first stage. The assumption \(\rho_m\to0\) but \(\rho_m\succ m^{-1/2}(\log m)^{1/2}\) is used to accommodate the difference between \(\boldsymbol{\theta}\) and its substitute \(\hat{\boldsymbol{\theta}}_m\), which satisfies \({\mathsf{P}}_{\boldsymbol{\theta}}\{\sqrt{m} \| \hat{\boldsymbol{\theta}}_m-\boldsymbol{\theta} \| \geqslant r(\log m)^{1/2}\} = o(m^{-1})\) if \(m\) is sufficiently large, by standard exponential bounds involving moment generating functions.

As noted by [72], the adaptive test in Section 6.2.2 can be regarded as a mid-course amendment of an adaptive test of \(H_0:u(\boldsymbol{\theta})\leqslant u_0\) versus \(H_1:u(\boldsymbol{\theta})\geqslant u_1\), with a maximum sample size of \(M\), to that of \(H_0\) versus \(H_2:u(\boldsymbol{\theta})\geqslant u_2\), with a maximum sample size of \(\widetilde{M}\). Whereas 82 provides an asymptotic lower bound for tests of \(H_0\) versus \(H_1\), any test of \(H_0\) versus \(H_2\) with error probabilities not exceeding \(\alpha\) and \(\tilde{\alpha}\) and taking at least \(m\) and at most \(\widetilde{M}\) observations likewise satisfies \[{\mathsf{E}}_{\boldsymbol{\theta}}(T)\geqslant m\vee\left\{\widetilde{M}\wedge\frac{[1+o(1)]|\log\alpha|}{I_0(\boldsymbol{\theta})\vee I_2(\boldsymbol{\theta})}\right\} \label{4466}\tag{83}\] as \(\alpha+\tilde{\alpha}\to0\) such that \(\log\alpha\sim\log\tilde{\alpha}\). Note that \(\Theta_1=\{\boldsymbol{\theta}:u(\boldsymbol{\theta})\geqslant u_1\} \subset\Theta_2=\{\boldsymbol{\theta}:u(\boldsymbol{\theta})\geqslant u_2\}\) and therefore \(I_2(\boldsymbol{\theta})\leqslant I_1(\boldsymbol{\theta})\). The four-stage test in Section 6.2.2, with \(M'=\widetilde{M}\), attempts to attain the asymptotic lower bound in 82 prior to the third stage and the asymptotic lower bound in 83 afterwards. It replaces \(I_1(\boldsymbol{\theta})\) in 82 , which corresponds to early stopping for futility, by \(I_2(\boldsymbol{\theta})\) that corresponds to rejection of \(H_2\) (instead of \(H_1\)) in favor of \(H_0\). Thus, the second-stage sample size \(n_2\) corresponds to the lower bound in 82 with \(\boldsymbol{\theta}\) replaced by \(\hat{\boldsymbol{\theta}}_m\) and \(I_1\) replaced by \(I_2\), while the third-stage sample size corresponds to that in 83 with \(\boldsymbol{\theta}\) replaced by \(\hat{\boldsymbol{\theta}}_{n_2}\). The arguments used to prove the asymptotic optimality of the three-stage test in Theorem 19 are then modified to prove the following.

Acknowledgements↩︎

JB: Lai was my postdoc advisor at Stanford in the “mid-aughts,” when I first got to witness his genius and work pace (“The speed of Lai,” as his Stanford colleague David Rogosa described it), usually in the Stanford Math Library in the middle of the night. I am grateful for Lai’s generosity in our research collaborations, his guidance of my young career, his contagious positivity, and the many laughs we shared along the way.

CDF: I have greatly enjoyed collaborating with Professor Tze-Liang Lai on Markov renewal theory. I am also deeply grateful to Prof. Lai for many fruitful conversations between 1990 and 2023, which played a significant role in my research, particularly in the areas of importance sampling and the multi-armed bandit problem.

AT: I am deeply grateful to Tze Lai for the many insightful conversations we shared between 1993 and 2023, which played a significant role in shaping my work. Tze’s work meaningfully influenced my research, from 1981 on.

HX: I am truly indebted to my advisor and collaborator, Professor Tze Leung Lai, for his guidance and inspiring mentorship. Our many discussions and collaboration over the years has profoundly shaped my research and the mentoring of students.

References↩︎