Papers

ergodicity is an open-source Python library for computational work on stochastic dynamics, with particular emphasis on non-ergodicity, time-average behavior, heavy-tailed processes, and decision making under uncertainty. The package brings together three layers that are often split across ad hoc scripts: process definitions and simulators, analysis and fitting tools, and agent-based experimentation. This article documents the implemented software rather than presenting new stochastic theory. We describe the package architecture, the supported process families, the analysis workflow, and the practical boundaries of the current implementation. We also provide fully reproducible examples covering heavy-tailed ensemble spread, multiplicative Levy growth diagnostics, adaptive memory mean reversion, preasymptotic fluctuation analysis, and partial stochastic differential equation simulation. The package is positioned as an integration layer on top of the scientific Python stack, reducing the amount of glue code required to move from process specification to diagnostics and comparative experiments.

Causal discovery, the problem of inferring the direction of causality, is generally ill-posed. We use the language of structural causal models (SCM) to show that assuming that the causal relations are acyclic and invariant across multiple environments (e.g., the way minimum wage affects employment rate is stable across different geographical regions), only two auxiliary environments are sufficient to infer the causal graph for arbitrary nonlinear mechanisms. Moreover, we demonstrate that this implies identifiability of the SCM functional mechanisms: as a corollary, we show that two auxiliary environments are sufficient to guarantee correct counterfactual inference. We empirically support our theoretical results on synthetic data.

Preprocessing screening is often the most expensive part of a near-infrared spectroscopy calibration workflow. It works because smoothing, derivatives, detrending and related filters change the spectral directions seen by partial least squares (PLS) or Ridge regression, but a full external search repeatedly refits nearly the same linear model. This paper studies the case where that search can be collapsed into one calibration step. For a strict linear preprocessing operator A acting on row spectra as XA^T, the transformed PLS cross-covariance satisfies (XA^T)^T Y = A X^T Y, and Ridge regression depends on the operator-induced kernel X A^T A X^T. These identities let a finite operator bank be screened inside the model while retaining original-wavelength coefficients, and the same identity extends to cheaply evaluated linear operator chains. Sample-adaptive or fitted corrections such as SNV, MSC, EMSC and ASLS are not strict linear; we prove the boundary and keep them as fold-local branches. The cohort has 61 regression and 17 classification rows, with a strict paired regression denominator of N=32 for the eight paper variants. There, AOM-PLS reaches median RMSEP ratios of 0.991/0.990 (simple) and 0.985/1.002 (best) against PLS-default/PLS-HPO, and AOM-Ridge reaches 0.974/0.984 (simple) and 0.918/0.966 (best) against Ridge-default/Ridge-HPO. The operator-adaptive classifier AOM-PLS-DA improves balanced accuracy by a median 0.159 on N=13 datasets (12/13 wins). The practical result is the runtime gap: PLS-HPO takes a median 710.81 s per run, whereas AOM-PLS takes 1.18-1.63 s – 436 to 602 times less PLS fitting time. Linear operator-adaptive calibration thus gives prediction quality comparable to exhaustive preprocessing screening, with orders-of-magnitude less fitting time for PLS.

Causal discovery methods aim to infer causal direction from observational data. Functional causal discovery approaches use structural asymmetries to identify causal directionality but rely on strong modeling assumptions and provide limited tools for uncertainty quantification. We introduce Causal Discovery via Statistical Power (CDSP), a statistical inference framework that connects causal direction estimation with statistical power and enables uncertainty quantification. Considering the foundational setting of bivariate observational data, we show how quantities analogous to statistical power and effect size can be used in causal discovery to determine when data contain sufficient information to favor one direction over the other. We introduce the effect-size asymmetry assumption that characterizes when the probability of correctly detecting the causal direction (i.e., the power of causal discovery) exceeds that of incorrectly favoring the reverse direction. We show that the effect-size asymmetry assumption can be used for causal direction estimation with uncertainty quantification. Simulations show that CDSP direction estimation is robust to mild and moderate model misspecifications. Real data analyses on 100 cause-effect benchmark pairs further demonstrate that CDSP reduces false discovery rates by approximately 18% relative to a commonly used existing method.

We investigate a geometric and distributional reinterpretation of Chatterjee's $\xi$-coefficient, which measures functional dependence between a response variable $Y$ and a predictor vector $\mathbf{X}$. For this purpose, we analyze the Markov product $(Y,Y')$, where $Y'$ is a copy of $Y$ that is conditionally independent of $Y$ given $\mathbf{X}$. Based on this construction, we introduce and study two dependence functions, denoted by $\phi_{(Y,\mathbf{X})}$ and $\kappa_{(Y,\mathbf{X})}$. The proposed framework provides a geometric interpretation of the Markov product and extends Chatterjee's correlation coefficient to a richer and more interpretable object for the analysis of directed stochastic dependence. In particular, rather than only measuring how well $Y$ can be represented as a function of $\mathbf{X}$, the proposed dependence functions additionally quantify how strongly the corresponding Markov product is concentrated near the diagonal.

The increasing availability of experimental data has intensified interest in calibrating stochastic models, raising fundamental questions about parameter identifiability. Structural identifiability determines whether parameters can be uniquely recovered from idealised, noise-free data, a prerequisite to allow for parameter estimation in real-world scenarios. However, existing methods to assess structural identifiability are not generally applicable to stochastic processes. We develop a methodology to analyse structural identifiability for a class of stochastic processes. We investigate how structural identifiability depends on the type of available data, distinguishing between single-particle trajectories and total particle density measurements. For trajectory data, we use the particle-based model description that explicitly represents single-particle dynamics. For population-level data, we derive a partial differential equation model representation, that describes the evolution of total particle density, and apply a differential algebra approach, common to ordinary differential equation analysis. We further introduce a method to study information arising from the initial condition, based on using the characteristic equations to construct a Taylor expansion of the particle density evolution. We apply our methodology to an example model and show that it is structurally identifiable from single-particle trajectory data but not from total particle density data, demonstrating that parameter identifiability depends on the type of data available.

Calibration is commonly evaluated by comparing model confidence with its empirical correctness, implicitly treating reliability as a function of the confidence score alone. However, this view can hide substantial structure: models may be systematically overconfident on some kinds of inputs and underconfident on others, causing global reliability diagnostics to obscure localised calibration failures. To address this, we formulate the problem of discovering hidden miscalibration regimes without assuming access to predefined data slices. We define the corresponding miscalibration field and propose a diagnostic framework for estimating it. Our approach learns a calibration-aware representation of the input space and estimates signed local miscalibration by kernel smoothing in the learned geometry. Across four real-world LLM benchmarks and twelve LLMs, we find that input-dependent calibration heterogeneity is prevalent. We further show that the discovered fields are actionable: they support local confidence correction and reduce calibration error in systematically miscalibrated regions where confidence-based methods such as isotonic regression and temperature scaling are less effective.

Diffusion models are often trained in low-dimensional latent spaces, which are then reused for related but shifted datasets. In this work, we study when such latent reuse remains reliable under distribution shift. We consider a source-target setting in which both datasets are approximately low-dimensional but may lie near different subspaces. We show that freezing and reusing a source latent space induces a target-domain score error governed by two quantities: the principal-angle misalignment between the source and target subspaces, and the target ambient noise amplified by the diffusion time scale. Motivated by these limits, we further study mixed source-target training and characterize how the required shared latent dimension depends on the relative geometry of the two distributions. Our results provide theoretical guidance on when latent reuse is reliable and when learning a shared representation may be necessary.

Continuous intraday electricity markets play an increasingly important role in short-term trading and balancing, yet decision-making under rapidly evolving price dynamics remains challenging. This paper proposes a comprehensive framework for ensemble forecasting of intraday electricity price trajectories and their translation into adaptive trading decisions. Building on a corrected Support Vector Regression model, the approach extends point predictions to probabilistic trajectory forecasts by introducing scenario generation based on forecast errors of fundamental variables and proposing a novel Support Vector Sorting procedure for the efficient selection of representative scenarios. The framework is evaluated using transaction level data from the German intraday continuous market. Empirical results show improvements over benchmark methods in both statistical and economic terms. Fundamental scenarios enhance median trajectory accuracy but produce more concentrated predictive distributions, while historical simulation with scenario selection better captures tail risk. From an economic perspective, ensemble-based forecasts outperform naive benchmarks across most of the trading strategies. Dynamic updating through scenario reweighting further improves profitability with limited impact on downside risk. Overall, the results demonstrate that combining kernel-based learning with scenario driven uncertainty and adaptive updating provides a flexible and effective approach for forecasting and trading in continuous electricity markets.

A median-radius framework for assessing centrality in multivariate data using median distances is proposed. Based on the proposed framework, a scale invariant measure of radial dispersion is defined and used to establish a depth function that is robust to outliers and independent of covariance structure. The depth function does not depend on moment assumptions and naturally adapts to skewness, multimodality, and heavy-tailed distributions, which make it effective for high-dimensional data structures. We demonstrate fundamental characteristics of the underlying functionals such as subgradient and convexity. The subgradients provide additional insight and encode the imbalance in directional contributions of the data. This suggests a new approach to detect skewness and structural asymmetry through a purely radial construction. Empirical studies demonstrate that the method agrees with classical approaches under symmetry while providing a more flexible and informative characterization in complex settings.

Asynchronous stochastic gradient descent (ASGD) is a standard way to exploit heterogeneous compute resources in distributed learning: instead of forcing fast workers to wait for slow ones, the server updates the model whenever a gradient arrives. Vanilla ASGD applies each arriving gradient with the same weight. When local data distributions are heterogeneous, this becomes problematic: faster workers contribute more updates, and we show theoretically that the method is biased toward a frequency-weighted average of the local objectives rather than the desired global objective. Existing remedies typically move away from the simple ASGD template by introducing gathering phases, buffering, or extra memory. We show that this is unnecessary. Keeping the standard ASGD mechanism, we recover the correct objective by rescaling worker-specific stepsizes in proportion to their computation times, so that each worker contributes the same aggregate learning rate over a cycle. In the non-convex setting, under smoothness and bounded heterogeneity assumptions, we prove that the resulting method, Rescaled ASGD, converges to stationary points of the correct global objective in the fixed-computation model. Its time complexity matches the known lower bound in the leading term, while the effects of staleness and data heterogeneity appear only in lower-order terms. Experiments confirm that the method converges to the correct objective and is competitive with state-of-the-art baselines.

Selection bias is pervasive in observational studies. For example, large scale biobanks data can exhibit “healthy volunteer bias” when respondents are healthier and of higher socio-economic status than the population they are meant to represent. Recovering causal effects from such sub-population is an important problem in causal inference, as estimating average treatment effects (ATE) from selected populations can result in a severely biased estimate of the ATE from the whole population. In this paper, we investigate the identifiability of the ATE under selection bias. We provide necessary and sufficient conditions for ATE identifiability, leveraging weak assumptions on probability classes to characterize propensity score and selection probability. Compared to previous works, our results extend existing graphical identifiability criteria and offer a more comprehensive understanding of causal effect identification with strictly weaker conditions in the presence of selection bias.

Many pre-trained models (PTMs) are available in modern applications. Because different PTMs are often trained on different datasets, their performances can vary substantially for different new tasks, and the ranking of the candidates may depend heavily on the input. Motivated by this, we propose a localized model averaging method with weights modeled as functions of the covariates, making it substantially more versatile than existing model averaging methods. This formulation allows the model averaging procedure to adaptively capture the varying relative advantages of different PTMs across heterogeneous contexts. Specifically, we learn flexible local weights under a general loss framework that accommodates a broad class of prediction tasks. We further establish the asymptotic optimality of the proposed method for both in-sample and out-of-sample risks, as well as the consistency of the estimated weights. Extensive numerical experiments further demonstrate the effectiveness of the proposed method.

We propose a data augmentation method for offline reinforcement learning, motivated by active positioning problems. Particularly, our approach enables the training of off-policy models from a limited number of suboptimal trajectories. We introduce a trajectory-based augmentation technique that exploits task structure and the geometric relationship between rewards, value functions, and mathematical properties of logging policies. During data collection, our augmentation supports suboptimal logging policies, leading to higher data quality and improved offline reinforcement learning performance. We provide theoretical justification for these strategies and validate them empirically across positioning tasks of varying dimensionality and under partial observability.

Bayesian inference for models with recursively defined likelihoods is computationally demanding, limiting scalability to large datasets. We propose a stabilised weighted subsampling methodology for accelerated Bayesian inference based on an unbiased estimator of the log-likelihood. By assigning higher sampling probabilities to early observations, the method reduces the effective depth of recursive likelihood evaluations and hence computational cost. This creates a challenging variance-cost trade-off: sampling probabilities that decay too slowly yield limited computational savings, while overly aggressive decay can substantially inflate estimator variance. We develop a stabilisation framework, supported by theory, that controls this trade-off and leads to a practical, theory-informed tuning strategy. We also derive an unbiased subsampling estimator of the log-likelihood gradient, enabling gradient-based Bayesian inference. The methodology can be embedded within a range of Bayesian computational algorithms. We illustrate its use in variational Bayes and subsampling Markov chain Monte Carlo for conditional volatility models, including leverage effects. Empirical results show substantial computational speed-ups relative to full-data methods while maintaining inferential accuracy. We also compare with recent stochastic gradient MCMC and divide-and-conquer MCMC methods for temporally dependent data, observing favourable empirical performance.

Applied researchers in biomedicine and related fields are often interested in estimating the causal effect of a treatment or intervention. Although randomized clinical trials are considered the gold standard for establishing causal effects, they are not always feasible, and real-world data may represent the only available source of evidence. In such settings, causal effects must be estimated using statistical methods applied to observational data. Over the last few decades, modern causal inference methods based on the potential outcomes framework have emerged as useful tools in this field. However, many such techniques exist, and their performance depends on factors such as sample size, the proportion of treated patients, the proportion of patients experiencing the outcome, the magnitude of the treatment effect, the target estimand, and potential violations of the fundamental assumptions of causal inference. Given the wide range of available methods, selecting an appropriate approach can be challenging for applied researchers. This study uses a large-scale simulation experiment to address this issue and provide researchers with a guide in the form of a handbook for a binary treatment and a binary outcome. Particularly, we test four popular statistical techniques: propensity score matching (full matching), inverse of the probability weighting, G-computation, and targeted maximum likelihood estimation. The proposed handbook is applied to two real-world datasets to assess its practical utility: one comprising vulnerable patients with mild COVID-19 (n=534 patients and more than 50% treated), and another of patients undergoing colorectal surgery (n=3635 patients and about 20% treated).

Generative models are often conditioned on a small set of examples via cross-attention. Under the Gaussian optimal-transport path, we show that the exact velocity field induced by a finite support set is a Nadaraya–Watson kernel smoother whose bandwidth decreases with flow time, from broad averaging at early steps to nearest-neighbor at late steps. A single Gaussian-kernel attention head exactly computes this field, connecting cross-attention conditioning to classical kernel theory. The theory predicts three failure regimes: nearest-neighbor collapse of the kernel at high dimension, mismatch between the isotropic kernel and the data geometry, and insufficient support for nonparametric estimation. Experiments on Gaussian mixtures, spherical shells, and DINOv2 ImageNet features confirm that learned conditioning improves in precisely these regimes, and that IP-Adapter's cross-attention implements approximate NW smoothing in practice.

This note addresses a key limitation of the Folding Test of Unimodality (FTU). In specific univariate mixture settings, the folding-based criterion can systematically fail, misclassifying clearly multimodal distributions as unimodal. We fully characterize these failures for Dirac mixtures and extend the analysis to Gaussian mixtures. We then introduce a double-folding procedure that captures complementary information, leading to a new test, the Double Folding Test of Unimodality. It resolves the FTU failures and improves multimodality detection power in simulations.

Most exploration algorithms search broadly until uncertainty is resolved. When the action space is too large to resolve within budget, practitioners default to $\varepsilon$-greedy, which bounds disruption but spends its override blindly. We introduce Delight-gated exploration (DE), a host–override rule that spends exploratory actions only when their prospective delight (expected improvement times surprisal) exceeds a gate price. This practical heuristic recovers a classical result: Pandora's reservation-value rule for costly search, with surprisal setting the effective inspection cost. Resolved arms exit the gate, fresh arms shut off above a prior-determined threshold, and selected linear-bandit overrides consume finite information budget. Across Bernoulli bandits, linear bandits, and tabular MDPs, the same hyperparameters transfer without retuning, and DE shows much weaker regret growth than Thompson Sampling and $\varepsilon$-greedy in the tested unresolved regimes. Delight improves acting for the same reason it improves learning: it prices scarce resources by the product of upside and surprisal.

Recent advancements in large language models demonstrate that injecting perturbations can substantially enhance extrapolation performance. However, current approaches often rely on discrete perturbations with fixed designs, which limits their flexibility. In this work, we propose a framework where token prefixes are perturbed by a learnable transformation of a continuous latent vector within an embedding space. To overcome the challenge of an intractable marginal likelihood, we derive unbiased estimating equations for model parameters and optimize them via stochastic gradient descent. We establish the statistical properties of the resulting estimator in over-parameterized regimes. Empirical evaluations on both synthetic and real-world datasets demonstrate that our proposal yields significant gains in out-of-domain settings over a range of state-of-the-art baseline methods.

We revisit Byzantine robust distributed estimation for high-dimensional sparse linear models. By combining local $\ell_1$-regularized robust estimation with robust aggregation at the server, the framework applies to pseudo-Huber regression, quantile regression, and sparse SVM. We show that the resulting estimators yield non-asymptotic guarantees and attain near-optimal statistical rates under mild conditions, while remaining communication-efficient. Simulations confirm strong robustness in estimation, support recovery and classification accuracy under various Byzantine attacks.

Physics-Informed Neural Networks (PINNs) and their variational counterparts (VPINNs) are neural networks that incorporate physical laws, making them useful for scientific problems. Existing generalization analyses for PINNs and VPINNs remain limited, often requiring restrictive assumptions such as stability conditions or linear ellipticity. In this paper, we derive generalization bounds for neural networks that involve differentiation with respect to input variables, covering PINNs and VPINNs under a unified framework. We apply Taylor expansion to represent nonlinear differential operators as linear operators on a high-dimensional space, enabling the use of Koopman-based analysis and showing that high-rank networks can generalize well even in settings involving differential operators. We also show that the nonlinearity of the differential operator exponentially enlarges the bound, highlighting its significant impact on generalization.

We study multiple change point localization under bandit feedback. An unknown piecewise-constant function on a compact interval can be queried sequentially at adaptively chosen inputs, and each query returns a noisy evaluation of the function. The goal is to identify a prescribed number of discontinuities, known as change points, within a target precision $\eta$ and confidence level $1-\delta$, while using as few samples as possible. We propose an adaptive algorithm that first detects intervals likely to contain change points and then refines their locations to precision $\eta$. We establish non-asymptotic upper bounds on its sample budget, together with corresponding lower bounds. Prior work shows that jump magnitudes alone determine the asymptotic sample complexity as $\delta\to 0$. We reveal that this picture is incomplete beyond this regime. We demonstrate, both empirically and theoretically, that for general $\delta$ and $\eta$, the complexity is jointly governed by the jumps and the relative positions of the change points.

This paper investigates the predictive performance of high-dimensional multimodel prediction, where the number of regressors is comparable to the sample size. Leveraging tools from random matrix theory, we derive the exact limiting out-of-sample risk under a nested model setting and comprehensively characterize the risk landscape. This limiting risk helps to reveal two phenomena: simple weighting inherits the double descent trajectory and its associated variance explosion near the interpolation boundary; strategic weighting triggers an ensemble emergence that suppresses the localized risk surge and yields a globally flat risk surface. Building on this limiting risk, we also propose the Large Model Averaging (LaMA) method, in which we consider the discrepancy between in-sample and out-of-sample risks in the high-dimensional regime. Numerical studies and real data applications confirm that LaMA achieves superior predictive accuracy in high-dimensional environments.

Large language models (LLMs) are increasingly deployed in settings where the available context is incomplete or degraded. We argue that an LLM generating answers under incomplete context can be viewed as an implicit imputer, and evaluated against a criterion from the multiple imputation (MI) literature: uncertainty should scale with the amount of missing information. We assess this criterion on SQuAD, using a controlled framework in which context availability is varied across five levels. We evaluate two answer-level uncertainty measures that can be estimated from repeated sampling: sampling-based confidence (empirical mode frequency) and response entropy. Confidence fails to reflect increasing missingness: it remains high even as accuracy collapses. Entropy, by contrast, increases with context removal, consistent with the MI analogy, and explains substantially more variance in accuracy than confidence across all evidence levels (quadratic $R^2$ gap up to 0.057). We further introduce a black-box diagnostic $\rho_R(\alpha)$ that estimates the proportion of baseline uncertainty resolved by context level $\alpha$, requiring only repeated sampling with and without context. These results suggest that entropy is a more responsive black-box uncertainty measure than confidence under incomplete context.

This work proposes $\chi^2$-type test statistics to assess different hypotheses on the local structure of an observed marked point pattern. The test statistics is based on the local inhomogeneous extension of the mark-weighted $K$-function to investigate local behaviour of the marked point pattern. The summary statistic captures interactions between marks and locations by assessing local contributions to global deviations from independence or homogeneity. The methodology proves to be effective in identifying both global and localised departures from the null hypotheses, even in scenarios with subtle mark structures or small sample sizes. Real-world environmental applications to forestry and earthquake data demonstrate the utility of the proposed framework for detecting spatially dependent marked structures in the patterns.

A likelihood-free transport filtering method is proposed based on the couplings between state and observation variables. By exploiting a block-triangular structure in the transport map, the analysis step of filtering is reformulated as the minimization of the maximum mean discrepancy (MMD) between the true joint measure and its transport-based approximation. To circumvent the non-convexity in the MMD optimization, we introduce a training-free transport filter method via gradient flows, which leads to an analytic computation for the transport map that implies the steepest descent direction of the MMD. The proposed approach accurately approximates non-Gaussian filtering posteriors and avoids particle collapse. We provide a convergence analysis for the expectation of the MMD between the approximated posterior and the truth posterior. Finally, we extend the method to high-dimensional problems through domain localization. Numerical examples demonstrate the superior performance of our approach over conventional filtering methods in nonlinear, non-Gaussian scenarios.

Michaelis–Menten analysis is often conducted by nonlinear least squares under a constant-variance assumption, even though enzyme-kinetic data frequently display concentration-dependent heteroscedasticity and often include repeated or clustered measurements. We develop a variance-aware procedure for Michaelis–Menten estimation and inference that is motivated by conditional moment restrictions and implemented through simple conditionally Gaussian working models. For single curves, the method reduces to one-dimensional root finding for $K_m$ followed by closed-form plug-in updates for $V_{\max}$ and a variance scale parameter; the same score logic yields a cluster-level extension through a random-effect-induced working covariance. In simulation, modeling heteroscedasticity improved variance recovery and interval efficiency relative to homoscedastic nonlinear least squares, while cluster-aware semiparametric and NLME fits restored fixed-effect coverage far more effectively than pooled analyses that ignored clustering. In self-driving laboratory and soil exoenzyme data, heteroscedastic models achieved lower information criteria than homoscedastic nonlinear least squares, with the square-root variance function giving the most stable empirical fit among the prespecified working models. We implement the workflow in the companion inferMM package for single-curve, grouped, and clustered Michaelis–Menten analysis. These results show that simple variance-function and covariance modeling can stabilize original-scale Michaelis–Menten inference when variability changes with substrate concentration or measurements are clustered.

Modern Bayesian optimization and adaptive sampling methods increasingly rely on nonlinear parametric models, yet theoretical guarantees for such models under adaptive data collection remain limited. Existing analyses largely focus on Gaussian processes, kernel machines, linear models, or linearized neural approximations, leaving a gap between theory and the nonlinear models used in practice. We develop a kernel based framework for analyzing regularized nonlinear parametric models trained on adaptively collected data. Our approach uses kernels over the parameter space to induce reproducing kernel Hilbert space structures over the corresponding model class, yielding confidence bounds for models trained with broad classes of regularized convex losses. We show how these bounds can support convergence guarantees for nonlinear acquisition and surrogate models, including randomized regularized policies that select points by maximizing a trained random model. These results provide a unified route to analyzing nonlinear parametric models in Bayesian optimization and related adaptive optimization settings.

Bell's theorem states that no description of a Bell experiment can be simultaneously local, realistic in the sense of counterfactual definiteness, and free of conspiracy between settings and hidden state. The recent generation of experiments has confirmed the predicted violation of the CHSH inequality, so one of the assumptions must be abandoned. Which one, and how one reconstructs a coherent worldview after doing so, is a question on which many authors disagree. This paper is written by three such authors. All three reject both counterfactual definiteness and conspiratorial violation of statistical independence of setting choices and state. After a joint exposition of the classical half of Bell's theorem in the language of Pearl-style causal graphs, a joint summary of the loophole-free experiments, and a joint survey of the recent literature, each author states where they have presently arrived. Gill accepts irreducible and non-local quantum randomness and finds the choice between locality and realism a false dichotomy. In his later works, Bell derives counterfactual definiteness from classical local causality, and that is what has to go. The metaphysical concepts "realism", "locality", "causality" need to be reconsidered. Helland reconstructs the Hilbert-space formalism from a theory of accessible variables, and from this theory he concludes that every observer must be limited in a specific sense. Jongejan proposes a geometric hidden-variable construction in which the degree of violation of the CHSH inequality depends on the number of dimensions of space, Tsirelson's bound corresponding to three dimensions. The authors conclude with a discussion.