Papers

Spherically embedded spatial data are spatially indexed observations whose values naturally reside on or can be equivalently mapped to the unit sphere. Such data are increasingly ubiquitous in fields ranging from geochemistry to demography. However, analysing such data presents unique difficulties due to the intrinsic non-Euclidean nature of the sphere, and rigorous methodologies for statistical modelling, inference, and uncertainty quantification remain limited. This paper introduces a unified framework to address these three limitations for spherically embedded spatial data. We first propose a novel spherical spatial autoregressive model that leverages optimal transport geometry and then extend it to accommodate exogenous covariates. Second, for either scenario with or without covariates, we establish the asymptotic properties of the estimators and derive a distribution-free Wald test for spatial dependence, complemented by a bootstrap procedure to enhance finite-sample performance. Third, we contribute a novel approach to uncertainty quantification by developing a conformal prediction procedure specifically tailored to spherically embedded spatial data. The practical utility of these methodological advances is illustrated through extensive simulations and applications to Spanish geochemical compositions and Japanese age-at-death mortality distributions.

Prediction markets are often described as mechanisms that “aggregate information” into prices, yet the mapping from dispersed private information to observed market histories is typically noisy, endogenous, and shaped by heterogeneous and strategic participation. This paper formulates prediction markets as Bayesian inverse problems in which the unknown event outcome \(Y\in\{0,1\}\) is inferred from an observed history of market-implied probabilities and traded volumes. We introduce a mechanism-agnostic observation model in log-odds space in which price increments conditional on volume arise from a latent mixture of trader types. The resulting likelihood class encompasses informed and uninformed trading, heavy-tailed microstructure noise, and adversarial or manipulative flow, while requiring only price and volume as observables. Within this framework we define posterior uncertainty quantification for \(Y\), provide identifiability and well-posedness criteria in terms of Kullback–Leibler separation between outcome-conditional increment laws, and derive posterior concentration statements and finite-sample error bounds under general regularity assumptions. We further study stability of posterior odds to perturbations of the observed price–volume path and define realized and expected information gain via the posterior-vs-prior KL divergence and mutual information. The inverse-problem formulation yields explicit diagnostics for regimes in which market histories are informative and stable versus regimes in which inference is ill-posed due to type-composition confounding or outcome–nuisance symmetries. Extensive experiments on synthetic data validate our theoretical predictions regarding posterior concentration rates and identifiability thresholds.

In competing risks models, cumulative incidence functions are commonly compared to infer differences between groups. Many existing inference methods, however, struggle when these functions cross during the time frame of interest. To address this problem, we investigate a test statistic based on the area between cumulative incidence functions. As the corresponding limiting distribution depends on quantities that are typically unknown, we propose a wild bootstrap approach to obtain a feasible and asymptotically valid two-sample test. The finite sample performance of the proposed method, in comparison with existing methods, is examined in an extensive simulation study.

Motivated by renal imaging studies that combine renogram curves with pharmacokinetic and demographic covariates, we propose Hybrid partial least squares (Hybrid PLS) for simultaneous supervised dimension reduction and regression in the presence of cross-modality correlations. The proposed approach embeds multiple functional and scalar predictors into a unified hybrid Hilbert space and rigorously extends the nonlinear iterative PLS (NIPALS) algorithm. This theoretical development is complemented by a sample-level algorithm that incorporates roughness penalties to control smoothness. By exploiting the rank-one structure of the resulting optimization problem, the algorithm admits a computationally efficient closed-form solution that requires solving only linear systems at each iteration. We establish fundamental geometric properties of the proposed framework, including orthogonality of the latent scores and PLS directions. Extensive numerical studies on synthetic data, together with an application to a renal imaging study, validate these theoretical results and demonstrate the method's ability to recover predictive structure under intermodal multicollinearity, yielding parsimonious low-dimensional representations.

Weather conditions can drastically alter the state of crops and rangelands, and in turn, impact the incomes and food security of individuals worldwide. Satellite-based remote sensing offers an effective way to monitor vegetation and climate variables on regional and global scales. The annual peak Normalized Difference Vegetation Index (NDVI), derived from satellite observations, is closely associated with crop development, rangeland biomass, and vegetation growth. Although various machine learning methods have been developed to forecast NDVI over short time ranges, such as one-month-ahead predictions, long-term forecasting approaches, such as one-year-ahead predictions of vegetation conditions, are not yet available. To fill this gap, we develop a two-phase machine learning model to forecast the one-year-ahead peak NDVI over high-resolution grids, using the Four Corners region of the Southwestern United States as a testbed. In phase one, we identify informative climate attributes, including precipitation and maximum vapor pressure deficit, and develop the generalized parallel Gaussian process that captures the relationship between climate attributes and NDVI. In phase two, we forecast these climate attributes using historical data at least one year before the NDVI prediction month, which then serve as inputs to forecast the peak NDVI at each spatial grid. We developed open-source tools that outperform alternative methods for both gross NDVI and grid-based NDVI one-year forecasts, providing information that can help farmers and ranchers make actionable plans a year in advance.

Longitudinal brain imaging data facilitate the monitoring of structural and functional alterations in individual brains across time, offering essential understanding of dynamic neurobiological mechanisms. Such data improve sensitivity for detecting early biomarkers of disease progression and enhance the evaluation of intervention effects. While recent matrix-response regression models can relate static brain networks to external predictors, there remain few statistical methods for longitudinal brain networks, especially those derived from high-dimensional imaging data. We introduce a matrix-response generalized linear mixed model that accommodates longitudinal brain networks and identifies edges whose connectivity is influenced by external predictors. An efficient Monte Carlo Expectation-Maximization algorithm is developed for parameter estimation. Extensive simulations demonstrate effective identification of covariate-related network components and accurate parameter estimation. We further demonstrate the usage of the proposed method through applications to diffusion tensor imaging (DTI) and functional MRI (fMRI) datasets.

It is recognised that treatment-related clustering should be allowed for in the sample size and analyses of individually-randomised parallel-group trials that evaluate therapist-delivered interventions such as psychotherapy. Here, interventions are a treatment factor, but therapists are not. If the aim of a trial is to separate effects of therapists from those of interventions, we propose that interventions and therapists should be regarded as two potentially interacting treatment factors (one fixed, one random) with a factorial structure. We consider the specific design where each therapist delivers each intervention (crossed therapist-intervention design), and the resulting therapist-intervention combinations are randomised to patients. We adopt a classical Design of Experiments (DoE) approach to propose a family of orthogonal factorial designs and their associated data analyses, which allow for therapist learning and centre too. We set out the associated data analyses using ANOVA and regression and report the results of a small simulation study conducted to explore the performance of the proposed randomisation methods in estimating the intervention effect and its standard error, the between-therapist variance and the between-therapist variance in the intervention effect. We conclude that more purposeful trial design has the potential to lead to better evidence on a range of complex interventions and outline areas for further methodological research.

We study a general framework of distributional computational graphs: computational graphs whose inputs are probability distributions rather than point values. We analyze the discretization error that arises when these graphs are evaluated using finite approximations of continuous probability distributions. Such an approximation might be the result of representing a continuous real-valued distribution using a discrete representation or from constructing an empirical distribution from samples (or might be the output of another distributional computational graph). We establish non-asymptotic error bounds in terms of the Wasserstein-1 distance, without imposing structural assumptions on the computational graph.

Despite the popular of multimodal statistical models, there lacks rigorous statistical inference tools for inferring the significance of a single modality within a multimodal model, especially in high-dimensional models. For high-dimensional multimodal generalized linear models, we propose a novel entropy-based metric, called the expected relative entropy, to quantify the information gain of one modality in addition to all other modalities in the model. We propose a deviance-based statistic to estimate the expected relative entropy, prove that it is consistent and its asymptotic distribution can be approximated by a non-central chi-squared distribution. That enables the calculation of confidence intervals and p-values to assess the significance of the expected relative entropy for a given modality. We numerically evaluate the empirical performance of our proposed inference tool by simulations and apply it to a multimodal neuroimaging dataset to demonstrate its good performance on various high-dimensional multimodal generalized linear models.

Machine-learned interatomic potentials (MLIPs) are increasingly used to replace computationally demanding electronic-structure calculations to model matter at the atomic scale. The most commonly used model architectures are constrained to fulfill a number of physical laws exactly, from geometric symmetries to energy conservation. Evidence is mounting that relaxing some of these constraints can be beneficial to the efficiency and (somewhat surprisingly) accuracy of MLIPs, even though care should be taken to avoid qualitative failures associated with the breaking of physical symmetries. Given the recent trend of scaling up models to larger numbers of parameters and training samples, a very important question is how unconstrained MLIPs behave in this limit. Here we investigate this issue, showing that – when trained on large datasets – unconstrained models can be superior in accuracy and speed when compared to physically constrained models. We assess these models both in terms of benchmark accuracy and in terms of usability in practical scenarios, focusing on static simulation workflows such as geometry optimization and lattice dynamics. We conclude that accurate unconstrained models can be applied with confidence, especially since simple inference-time modifications can be used to recover observables that are consistent with the relevant physical symmetries.

Uncertainty estimation in machine learning has traditionally focused on the prediction stage, aiming to quantify confidence in model outputs while treating learned representations as deterministic and reliable by default. In this work, we challenge this implicit assumption and argue that reliability should be regarded as a first-class property of learned representations themselves. We propose a principled framework for reliable representation learning that explicitly models representation-level uncertainty and leverages structural constraints as inductive biases to regularize the space of feasible representations. Our approach introduces uncertainty-aware regularization directly in the representation space, encouraging representations that are not only predictive but also stable, well-calibrated, and robust to noise and structural perturbations. Structural constraints, such as sparsity, relational structure, or feature-group dependencies, are incorporated to define meaningful geometry and reduce spurious variability in learned representations, without assuming fully correct or noise-free structure. Importantly, the proposed framework is independent of specific model architectures and can be integrated with a wide range of representation learning methods.

Imbalanced classification often causes standard training procedures to prioritize the majority class and perform poorly on rare but important cases. A classic and widely used remedy is to augment the minority class with synthetic samples, but two basic questions remain under-resolved: when does synthetic augmentation actually help, and how many synthetic samples should be generated? We develop a unified statistical framework for synthetic augmentation in imbalanced learning, studying models trained on imbalanced data augmented with synthetic minority samples. Our theory shows that synthetic data is not always beneficial. In a "local symmetry" regime, imbalance is not the dominant source of error, so adding synthetic samples cannot improve learning rates and can even degrade performance by amplifying generator mismatch. When augmentation can help ("local asymmetry"), the optimal synthetic size depends on generator accuracy and on whether the generator's residual mismatch is directionally aligned with the intrinsic majority-minority shift. This structure can make the best synthetic size deviate from naive full balancing. Practically, we recommend Validation-Tuned Synthetic Size (VTSS): select the synthetic size by minimizing balanced validation loss over a range centered near the fully balanced baseline, while allowing meaningful departures. Extensive simulations and real data analysis further support our findings.

In this paper, we study the spherical cardioid distribution, a higher-dimensional and arbitrary-order generalization of the circular cardioid distribution. This distribution is rotationally symmetric and generates unimodal, multimodal, axial, and girdle-like densities. We identify various properties of the spherical cardioid that make it highly tractable: simple density evaluation, closedness under convolution, explicit expressions for vectorized moments, and efficient simulation. The moments of the spherical cardioid of order $k$ up to order $k-1$ coincide with those of the uniform distribution on the sphere, highlighting its closeness to the latter. We derive estimators by the method of moments and maximum likelihood, their asymptotic distributions, and their asymptotic relative efficiencies. We give the machinery for bootstrap goodness-of-fit tests based on the projected empirical cumulative distribution function approach, including the projected distribution and closed-form expressions for test statistics. An application to modeling the orbits of long-period comets shows the usefulness of the spherical cardioid distribution in real data analyses.

We introduce a new online approach for constructing proposal distributions in particle filters using a forward scheme. Our method progressively incorporates future observations to refine proposals. This is in contrast to backward-scheme algorithms that require access to the entire dataset, such as the iterated auxiliary particle filters (Guarniero et al., 2017, arXiv:1511.06286) and controlled sequential Monte Carlo (Heng et al., 2020, arXiv:1708.08396) which leverage all future observations through backward recursion. In comparison, our forward scheme achieves a gradual improvement of proposals that converges toward the proposal targeted by these backward methods. We show that backward approaches can be numerically unstable even in simple settings. Our forward method, however, offers significantly greater robustness with only a minor trade-off in performance, measured by the variance of the marginal likelihood estimator. Numerical experiments on both simulated and real data illustrate the enhanced stability of our forward approach.

Deep neural networks (DNNs) typically involve a large number of parameters and are trained to achieve zero or near-zero training error. Despite such interpolation, they often exhibit strong generalization performance on unseen data, a phenomenon that has motivated extensive theoretical investigations. Comforting results show that interpolation indeed may not affect the minimax rate of convergence under the squared error loss. In the mean time, DNNs are well known to be highly vulnerable to adversarial perturbations in future inputs. A natural question then arises: Can interpolation also escape from suboptimal performance under a future $X$-attack? In this paper, we investigate the adversarial robustness of interpolating estimators in a framework of nonparametric regression. A finding is that interpolating estimators must be suboptimal even under a subtle future $X$-attack, and achieving perfect fitting can substantially damage their robustness. An interesting phenomenon in the high interpolation regime, which we term the curse of simple size, is also revealed and discussed. Numerical experiments support our theoretical findings.

Change point tests for abrupt changes in the mean of functional data, i.e., random elements in infinite-dimensional Hilbert spaces, are either based on dimension reduction techniques, e.g., based on principal components, or directly based on a functional CUSUM (cumulative sum) statistic. The former have often been criticized as not being fully functional and losing too much information. On the other hand, unlike the latter, they take the covariance structure of the data into account by weighting the CUSUM statistics obtained after dimension reduction with the inverse covariance matrix. In this paper, as a middle ground between these two approaches, we propose an alternative statistic that includes the covariance structure with an offset parameter to produce a scale-invariant test procedure and to increase power when the change is not aligned with the first components. We obtain the asymptotic distribution under the null hypothesis for this new test statistic, allowing for time dependence of the data. Furthermore, we introduce versions of all three test statistics for gradual change situations, which have not been previously considered for functional data, and derive their limit distribution. Further results shed light on the asymptotic power behavior for all test statistics under various ground truths for the alternatives.

In constrained stochastic optimization, one expects that restricting the feasible set, provided it still contains the true parameter, should not increase the statistical risk of the corresponding projection estimator. We show that this intuition can fail, even in basic settings. We investigate this phenomenon in the Gaussian sequence model. Given a compact, convex set $\Theta \subseteq \mathbb{R}^d$, one observes \[ Y = \theta^\star + \sigma Z, \qquad Z \sim N(0, I_d), \] and seeks to estimate an unknown $\theta^\star \in \Theta$. Here, the maximum likelihood estimator over $\Theta$ coincides with the least squares estimator (LSE), given by the Euclidean projection of $Y$ onto $\Theta$. We construct an explicit example exhibiting risk reversal: for sufficiently large noise, there exist nested compact convex sets $\Theta_S \subsetneq \Theta_L$ and $\theta^\star \in \Theta_S$ such that the LSE constrained to $\Theta_S$ has strictly larger squared-error risk than the LSE constrained to $\Theta_L$. Moreover, we demonstrate that risk reversal can persist at the level of worst-case risk. Finally, we show that the phenomenon is not specific to Gaussian noise or squared-error risk, extending our results beyond both settings. We clarify this phenomenon by contrasting noise regimes. In the vanishing-noise limit, the risk is governed at first order by the statistical dimension of the tangent cone, and risk reversal cannot occur at the leading $\sigma^2$ scale. In the diverging-noise regime, the risk instead depends on the global geometry of the constraint sets, and the embedding of $\Theta_S$ within $\Theta_L$ can reverse the risk ordering. These results reveal a previously unrecognized failure mode of the LSE. They demonstrate that in sufficiently noisy settings, tightening a constraint can paradoxically degrade statistical performance.

We study maximum likelihood estimation for spatial generalized linear mixed models with Gaussian process approximations using a stochastic Newton-Raphson algorithm. We consider two Gaussian Process approximations in this context: spectral Gaussian process approximations and stochastic partial differential equations (SPDE). We refine the stochastic maximum likelihood algorithm and we propose a new stopping criterion for efficient termination to prevent long runs of sampling in the stationary post-convergence phase and a Monte Carlo estimator of fixed effect standard errors. We run a series of simulation comparisons of spatial statistical models alongside the popular Bayesian integrated nested Laplacian approximation method which incorporates SPDE. We show that HSGP provides nominal coverage of fixed and random effect parameters with smooth latent fields but performance degrades for rough fields. SPDE in a stochastic maximum likelihood framework maintains nominal coverage and matches or improves upon the performance of Bayesian integrated nested Laplacian approximation.

This paper investigates the minimax-optimality of Halpern fixed-point iterations for Lipschitz maps in general normed spaces. Starting from an a priori bound on the orbit of iterates, we derive non-asymptotic estimates for the fixed-point residuals. These bounds are tight, meaning that they are attained by a suitable Lipschitz map and an associated Halpern sequence. By minimizing these tight bounds we identify the minimax-optimal Halpern scheme. For contractions, the optimal iteration exhibits a transition from an initial Halpern phase to the classical Banach-Picard iteration and, as the Lipschitz constant approaches one, we recover the known convergence rate for nonexpansive maps. For expansive maps, the algorithm is purely Halpern with no Banach-Picard phase; moreover, on bounded domains, the residual estimates converge to the minimal displacement bound. Inspired by the minimax-optimal iteration, we design an adaptive scheme whose residuals are uniformly smaller than the minimax-optimal bounds, and can be significantly sharper in practice. Finally, we extend the analysis by introducing alternative bounds based on the distance to a fixed point, which allow us to handle mappings on unbounded domains; including the case of affine maps for which we also identify the minimax-optimal iteration.

We consider a general class of round-robin tournament models of equally strong players. In these models, each of the $n$ players competes against every other player exactly once. For each match between two players, the outcome is a value from a countable subset of the unit interval, and the scores of the two players in a match sum to one. The final score of each player is defined as the sum of the scores obtained in matches against all other players. We study the distribution of extreme scores, including the maximum, second maximum, and lower-order extremes. Since the exact distribution is computationally intractable even for small values of $n$, we derive asymptotic results as the number of players $n$ tends to infinity, including limiting distributions, and rates of convergence.

In prognostics and health management (PHM) of engineered systems, maintenance decisions are ideally informed by predictions of a system's remaining useful life (RUL) based on operational data. Model-based prognostics algorithms rely on a parametric model of the system degradation process. The model parameters are learned from real-time operational data collected on the system. However, there can be valuable information in data from similar systems or components, which is not typically utilized in PHM. In this contribution, we propose a hierarchical Bayesian modeling (HBM) framework for PHM that integrates both operational data and run-to-failure data from similar systems or components. The HBM framework utilizes hyperparameter distributions learned from data of similar systems or components as priors. It enables efficient updates of predictions as more information becomes available, allowing for increasingly accurate assessments of the degradation process and its associated variability. The effectiveness of the proposed framework is demonstrated through two experimental applications involving real-world data from crack growth and lithium battery degradation. Results show significant improvements in RUL prediction accuracy and demonstrate how the framework facilitates uncertainty management through predictive distributions.

Understanding historical datasets, such as the England and Wales infant mortality data, for local government districts can provide valuable insights into our changing society. Such analyses can prove challenging in practice, due to frequent changes in the boundaries of local government districts for which records are collected. One solution adopted in the literature to overcome such practical challenges is to pre-process data using areal interpolation to render the units consistent over the time period of focus. However, such methods are prone to errors. In this paper we introduce a novel changepoint method to detect instances where interpolation performs poorly. We demonstrate the utility of our method on original data, and also demonstrate how correcting interpolation errors can affect the clustering of the infant mortality curves.

We study two-sample equality testing in Gaussian graphical models. Classical likelihood ratio tests on decomposable graphs admit clique-wise factorizations, offering limited localization and unstable finite-sample behaviour. We propose node-level inference via a leave-one-out Bartlett-adjusted test on a fully connected graph. The resulting increments have standard chi-square null limits, enabling calibrated significance for single nodes and fixed-size subsets. Simulations confirm validity, and a case study shows practical utility.

The Mann-Whitney effect is an effect measure for the order of two sample-specific outcome variables. It has the interpretation of a probability and also a connection to the area under the ROC curve. In the literature it has been considered for both ordinal and right-censored time-to-event outcomes. For both cases, the present paper introduces a distribution-free regression model that relates the Mann-Whitney effect to a linear combination of covariates. To fit the model, we develop a pseudo-observation-based procedure yielding consistent and asymptotically normal coefficient estimates. In addition, we propose bootstrap-based hypothesis tests to infer the effects of the covariates on the Mann-Whitney effect. A simulation study on the small-sample behavior of the proposed method demonstrates that the novel hypothesis tests keep up with the z-test of a Cox regression model. The new methods are used to analyze progression-free survival in breast cancer patients enrolled for the randomized phase III SUCCESS-A trial.

We introduce the AlgebraicStatistics section of the OSCAR computer algebra system. We give an overview of its extensible design and highlight its features including serialization of data types for sharing results and creating databases, and state-of-the-art implicitization algorithms.

Functional graphical models have undergone extensive development during the recent years, leading to a variety models such as the functional Gaussian graphical model, the functional copula Gaussian graphical model, the functional Bayesian graphical model, the nonparametric functional additive graphical model, and the conditional functional graphical model. These models rely either on some parametric form of distributions on random functions, or on additive conditional independence, a criterion that is different from probabilistic conditional independence. In this paper we introduce a nonparametric functional graphical model based on functional sufficient dimension reduction. Our method not only relaxes the Gaussian or copula Gaussian assumptions, but also enhances estimation accuracy by avoiding the “curse of dimensionality”. Moreover, it retains the probabilistic conditional independence as the criterion to determine the absence of edges. By doing simulation study and analysis of the f-MRI dataset, we demonstrate the advantages of our method.

While Large Language Models (LLMs) show remarkable capabilities, their unreliability remains a critical barrier to deployment in high-stakes domains. This survey charts a functional evolution in addressing this challenge: the evolution of uncertainty from a passive diagnostic metric to an active control signal guiding real-time model behavior. We demonstrate how uncertainty is leveraged as an active control signal across three frontiers: in advanced reasoning to optimize computation and trigger self-correction; in autonomous agents to govern metacognitive decisions about tool use and information seeking; and in reinforcement learning to mitigate reward hacking and enable self-improvement via intrinsic rewards. By grounding these advancements in emerging theoretical frameworks like Bayesian methods and Conformal Prediction, we provide a unified perspective on this transformative trend. This survey provides a comprehensive overview, critical analysis, and practical design patterns, arguing that mastering the new trend of uncertainty is essential for building the next generation of scalable, reliable, and trustworthy AI.

Bayesian optimisation is a sample efficient method for finding a global optimum of expensive black-box objective functions. Historic datasets from related problems can be exploited to help improve performance of Bayesian optimisation by adapting transfer learning methods to various components of the Bayesian optimisation pipeline. In this study we perform an empirical analysis of various ensemble-based transfer learning Bayesian optimisation methods and pipeline components. We expand on previous work in the literature by contributing some specific pipeline components, and three new real-time transfer learning Bayesian optimisation benchmarks. In particular we propose to use a weighting strategy for ensemble surrogate model predictions based on regularised regression with weights constrained to be positive, and a related component for handling the case when transfer learning is not improving Bayesian optimisation performance. We find that in general, two components that help improve transfer learning Bayesian optimisation performance are warm start initialisation and constraining weights used with ensemble surrogate model to be positive.

In the study of time-dependent (i.e., temporal) networks, researchers often examine the evolution of communities, which are sets of densely connected sets of nodes that are connected sparsely to other nodes. An increasingly prominent approach to studying community structure in temporal networks is statistical inference. In the present paper, we study the performance of a class of statistical-inference methods for community detection in temporal networks. We represent temporal networks as multilayer networks, with each layer encoding a time step, and we illustrate that statistical-inference models that generate community assignments via either a uniform distribution on community assignments or discrete-time Markov processes are biased against generating communities with large or small numbers of nodes. In particular, we demonstrate that statistical-inference methods that use such generative models tend to poorly identify community structure in networks with large or small communities. To rectify this issue, we introduce a novel statistical model that generates the community assignments of the nodes in given layer (i.e., at a given time) using all of the community assignments in the previous layer. We prove results that guarantee that our approach greatly mitigates the bias against large and small communities, so using our generative model is beneficial for studying community structure in networks with large or small communities. Our code is available at https://github.com/tfaust0196/TemporalCommunityComparison.

Recently, Major League Baseball (MLB) limited pitchers to three pickoff attempts, creating a cat-and-mouse game between pitcher and runner. Each failed attempt adds pressure on the pitcher to avoid using another, and the runner can intensify this pressure by extending their leadoff toward the next base. We model this dynamic as a two-player zero-sum sequential game in which the runner first chooses a lead distance, and then the pitcher chooses whether to attempt a pickoff. We establish optimality characterizations for the game and present variants of value iteration and policy iteration to solve the game. Using high-resolution lead distance data from optical player tracking, we estimate generalized linear mixed-effects models for pickoff and stolen base outcome probabilities given lead distance, context, and player skill. We compute the game-theoretic equilibria under the two-player model, as well as the optimal runner policy under a simplified one-player Markov decision process (MDP) model. In the one-player setting, our results establish an actionable rule of thumb: the Two-Foot Rule, which recommends that a runner increase their lead by two feet after each pickoff attempt.