Current protocols to estimate the number, size, and location of cancerous lesions in the prostate using multiparametric magnetic resonance imaging (mpMRI) are highly dependent on reader experience and expertise. Automatic voxel-wise cancer classifiers do not directly provide estimates of number, location, and size of cancerous lesions that are clinically important. Existing spatial partitioning methods estimate linear or piecewise-linear boundaries separating regions of local stationarity in spatially registered data and are inadequate for the application of lesion detection. Frequentist segmentation and clustering methods often require pre-specification of the number of clusters and do not quantify uncertainty. Previously, we developed a novel Bayesian functional spatial partitioning method to estimate the boundary surrounding a single cancerous lesion using data derived from mpMRI. We propose a Bayesian functional spatial partitioning method for multiple lesion detection with an unknown number of lesions. Our method utilizes functional estimation to model the smooth boundary curves surrounding each cancerous lesion. In a Reversible Jump Markov Chain Monte Carlo (RJ-MCMC) framework, we develop novel jump steps to jointly estimate and quantify uncertainty in the number of lesions, their boundaries, and the spatial parameters in each lesion. Through simulation we show that our method is robust to the shape of the lesions, number of lesions, and region-specific spatial processes. We illustrate our method through the detection of prostate cancer lesions using MRI.

When analysing multiple time series that may be subject to changepoints, it is sometimes possible to specify a priori, by means of a graph, which pairs of time series are likely to be impacted by simultaneous changepoints. This article proposes an informative prior for changepoints which encodes the information contained in the graph, inducing a changepoint model for multiple time series that borrows strength across clusters of connected time series to detect weak signals for synchronous changepoints. The graphical model for changepoints is further extended to allow dependence between nearby but not necessarily synchronous changepoints across neighbouring time series in the graph. A novel reversible jump Markov chain Monte Carlo (MCMC) algorithm making use of auxiliary variables is proposed to sample from the graphical changepoint model. The merit of the proposed approach is demonstrated through a changepoint analysis of computer network authentication logs from Los Alamos National Laboratory (LANL), demonstrating an improvement at detecting weak signals for network intrusions across users linked by network connectivity, whilst limiting the number of false alerts.

Approximate Bayesian computing (ABC) and Bayesian Synthetic likelihood (BSL) are two popular families of methods to evaluate the posterior distribution when the likelihood function is not available or tractable. For existing variants of ABC and BSL, the focus is usually first put on the simulation algorithm, and after that the form of the resulting approximate posterior distribution comes as a consequence of the algorithm. In this paper we turn this around and firstly define a reasonable approximate posterior distribution by studying the distributional properties of the expected discrepancy, or more generally an expected evaluation, with respect to generated samples from the model. The resulting approximate posterior distribution will be on a simple and interpretable form compared to ABC and BSL.

Secondly a Markov chain Monte Carlo (MCMC) algorithm is developed to simulate from the resulting approximate posterior distribution. The algorithm was evaluated on a synthetic data example and on the Stepping Stone population genetics model, demonstrating that the proposed scheme has real world applicability. The algorithm demonstrates competitive results with the BSL and sequential Monte Carlo ABC algorithms, but is outperformed by the ABC MCMC.

Testing the (in)equality of variances is an important problem in many statistical applications. We develop default Bayes factor tests to assess the (in)equality of two or more population variances, as well as a test for whether a population variance equals a specific value. The resulting test can be used to check assumptions for commonly used procedures such as the t-test or ANOVA, or test substantive hypotheses concerning variances directly. We show that our Bayes factor fulfills a number of desiderata. Researchers may have directed hypotheses such as ${\mathit{\sigma }}_1^2>{\mathit{\sigma }}_2^2$, they may want to extend $H_0$ to have a null-region, or wish to combine hypotheses about equality with hypotheses about inequality, for example ${\mathit{\sigma }}_1^2={\mathit{\sigma }}_2^2>({\mathit{\sigma }}_3^2,{\mathit{\sigma }}_4^2)$. We extend our Bayes factor test to allow for these deviations from our proposed default and illustrate it on a number of practical examples. Our procedure is implemented in the R package bfvartest.

The objective of disease mapping is to model data aggregated at the areal level. In some contexts, however, (e.g. residential histories, general practitioner catchment areas) when data is arising from a variety of sources, not necessarily at the same spatial scale, it is possible to specify spatial random effects, or covariate effects, at the areal level, by using a multiple membership principle (MM) (Petrof et al., 2020; Gramatica et al., 2021). A weighted average of conditional autoregressive (CAR) spatial random effects embeds spatial information for a spatially-misaligned outcome and estimate relative risk for both frameworks (areas and memberships). In this paper we investigate the theoretical underpinnings of these application of the multiple membership principle to the CAR prior, in particular with regard to parameterisation, properness and identifiability. We carry out simulations involving different numbers of memberships as compared to number of areas and assess impact of this on estimating parameters of interest. Both analytical and simulation study results show under which conditions parameters of interest are identifiable, so that we can offer actionable recommendations to researchers. Finally, we present the results of an application of the multiple membership model to diabetes prevalence data in South London, together with strategic implications for public health considerations.

This research is focused on a formal Bayesian method of recursive multi-step-ahead density prediction and its ex post evaluation. Our approach remains within the framework of the standard (classical or orthodox) Bayesian paradigm based on the Bayes factor and on the use of the likelihood-based update. We propose a new decomposition of the predictive Bayes factor into the product of partial Bayes factors, for both a finite number of consecutive k-step-ahead forecasts (where $k>1$) and the recursive updates of the posterior odds ratios based on updated data sets. The first factor in the decomposition is related to the relative k-step-ahead forecasting ability of models, while the second one measures the updating effect.

To illustrate the usefulness of the proposed measures, we apply the new decomposed predictive Bayes factors to compare the forecasting ability of models when the true data generating process (DGP) is known, using simulated data sets. Taking into account the effect of updating, the posterior odds ratios leads to the conclusion that the best model coincides with the true DGP. However, the highest k-step-ahead forecasting ability (considered alone) can be achieved by some other, less adequate models. Next, we investigate the predictive ability of different Vector Error Correction (VEC) models with conditional heteroscedasticity, combining three macroeconomic variables: unemployment, inflation and interest rates, separately for the US and Polish economies. The results show that the inference about the models’ predictive performance depends on the forecast horizon as well as on taking into account the updating effect.

Heavy-tailed continuous shrinkage priors, such as the horseshoe prior, are widely used for sparse estimation problems. However, there is limited work extending these priors to explicitly incorporate multivariate shrinkage for regressors with grouping structures. Of particular interest in this article, is regression coefficient estimation where pockets of high collinearity in the regressor space are contained within known regressor groupings. To assuage variance inflation due to multicollinearity we propose the group inverse-gamma gamma (GIGG) prior, a heavy-tailed prior that can trade-off between local and group shrinkage in a data adaptive fashion. A special case of the GIGG prior is the group horseshoe prior, whose shrinkage profile is dependent within-group such that the regression coefficients marginally have exact horseshoe regularization. We establish posterior consistency and posterior concentration results for regression coefficients in linear models and mean parameters in sparse normal means models. The full conditional distributions corresponding to GIGG regression can be derived in closed form, leading to straightforward posterior computation. We show that GIGG regression results in low mean-squared error across a wide range of correlation structures and within-group signal densities via simulation. We apply GIGG regression to data from the National Health and Nutrition Examination Survey for associating environmental exposures with liver functionality.

The efficiency of Hamiltonian Monte Carlo (HMC) can suffer when sampling a distribution with a wide range of length scales, because the small step sizes needed for stability in high-curvature regions are inefficient elsewhere. To address this we present a delayed rejection (DR) variant: if an initial HMC trajectory is rejected, we make one or more subsequent proposals each using a step size geometrically smaller than the last. To reduce the cost of DR approaches, we extend the standard delayed rejection to a probabilistic framework wherein we do not make multiple proposals at every rejection, but allow the probability of a retry to depend on the probability of accepting the previous proposal. We test the scheme in several sampling tasks, including statistical applications and multiscale model distributions such as Neal’s funnel. Delayed rejection enables sampling multiscale distributions for which standard approaches such as HMC fail to explore the tails, and improves performance five-fold over optimally-tuned HMC as measured by effective sample size per gradient evaluation. Even for simpler distributions, delayed rejection provides increased robustness to step size misspecification.

Studies with binary outcomes on a heterogeneous population are quite common. Typically, the heterogeneity is modelled through varying effect coefficients within some binary regression setting combined with a clustering procedure. Most of the existing methods assign statistical units to distinct and non-overlapping clusters. However, there are scenarios where units exhibit a more complex organization and the clusters can be thought as partially overlapping. In this case, the standard approach does not work. In this paper, we define a mixture of regression models that allows overlapping clusters. This approach involves an overlap function that maps the regression coefficients, either at the unit or response level, of the parent clusters into the coefficients of the multiple allocation clusters. In order to deal with this intrinsic heterogeneity, regression analyses have to be stratified for different groups of observations or clusters. We present a computationally efficient Monte Carlo Markov Chain (MCMC) scheme for the case of a mixture of probit regressions. A simulation study shows the overall performance of the method. We conclude with two illustrative examples of modelling voting behavior, involving United States (US) Supreme Court justices over a number of topics and members of the United Kingdom (UK) parliament over divisions related to Brexit. These applications provide insights on the usefulness of the method in real applications. The method described can be extended to the case of a generic mixture of multivariate generalized linear models under overlapping clusters.

We propose minimax optimal Bayesian two-sample tests for testing equality of high-dimensional mean vectors and covariance matrices between two populations. In many applications including genomics and medical imaging, it is natural to assume that only a few entries of two mean vectors or covariance matrices are different. Many existing tests that rely on aggregating the difference between empirical means or covariance matrices are not optimal or yield low power under such setups. Motivated by this, we develop Bayesian two-sample tests employing a divide-and-conquer idea, which is powerful especially when the differences between two populations are rare but large. The proposed two-sample tests manifest closed forms of Bayes factors and allow scalable computations even in high-dimensions. We prove that the proposed tests are consistent under relatively mild conditions compared to existing tests in the literature. Furthermore, the testable regions from the proposed tests turn out to be minimax optimal in terms of rates. Simulation studies show clear advantages of the proposed tests over other state-of-the-art methods in various scenarios. Our tests are also applied to the analysis of the gene expression data of two cancer data sets.

We consider the joint inference of regression coefficients and the inverse covariance matrix for covariates in high-dimensional probit regression, where the predictors are both relevant to the binary response and functionally related to one another. A hierarchical model with spike and slab priors over regression coefficients and the elements in the inverse covariance matrix is employed to simultaneously perform variable and graph selection. We establish joint selection consistency for both the variable and the underlying graph when the dimension of predictors is allowed to grow much larger than the sample size, which is the first theoretical result in the Bayesian literature. A scalable Gibbs sampler is derived that performs better in high-dimensional simulation studies compared with other state-of-art methods. We illustrate the practical impact and utilities of the proposed method via a functional MRI dataset, where both the regions of interest with altered functional activities and the underlying functional brain network are inferred and integrated together for stratifying disease risk.

We develop Bayesian model selection via fractional Bayes factors to simultaneously assess spatial dependence and select regressors in Gaussian hierarchical models with intrinsic conditional autoregressive (ICAR) spatial random effects. Selection of covariates and spatial model structure is difficult, as spatial confounding creates a tension between fixed and spatial random effects. Researchers have commonly performed selection separately for fixed and random effects in spatial hierarchical models. Simultaneous selection methods relieve the researcher from arbitrarily fixing one of these types of effects while selecting the other. Notably, Bayesian approaches to simultaneously select covariates and spatial model structure are limited. Our use of fractional Bayes factors allows for selection of fixed effects and spatial model structure under automatic reference priors for model parameters, which obviates the need to specify hyperparameters for priors. We also show the equivalence between two ICAR specifications and derive the minimal training size for the fractional Bayes factor applied to the ICAR model under the reference prior. We perform a simulation study to assess the performance of our approach and we compare results to the Deviance Information Criterion and Widely Applicable Information Criterion. We demonstrate that our fractional Bayes factor approach assigns low posterior model probability to spatial models when data is truly independent and reliably selects the correct covariate structure with highest probability within the model space. Finally, we demonstrate our Bayesian model selection approach with applications to county-level median household income in the contiguous United States and residential crime rates in the neighborhoods of Columbus, Ohio.

This expository paper discusses Bayesian decision analysis perspectives on problems of constrained forecasting. Foundational and pedagogic discussion contrasts decision analytic approaches with the traditional, but typically inappropriate, inferential approach. Illustrative examples include development of novel constrained point forecasting and entropic tilting methodology to explore consistency of a predictive distribution with an imposed or hypothesized constraint. Linear, aggregate constraints define illuminating examples that relate to broadly important problems involving aggregate and hierarchical constraints in commercial and economic forecasting. Discussion explores the impact of different loss functions, questions of how constrained forecasting is impacted by dependencies among outcomes being predicted, and promotes the broader use of decision analysis including routine evaluation of predictive distributions of loss under chosen forecasts/decisions. Extensions to more general constrained forecasting problems, connections with broader interests in forecast reconciliation and other considerations are noted.

This paper proposes a novel approach that combines dynamic linear models applied to graph data and variable selection through spike-and-slab priors. The new class of models, called Dynamic Graphical Variable Selection, is used to infer effective connectivity in certain brain regions allowing both connectivity weights and graphical structure to vary over time. One advantage of our method is that as the graphical structure is estimated inferentially, the computational cost is reduced. That way our methodology can accommodate high-dimensional data, such as large networks observed through long periods of time. We illustrate our methodology via numerical experiments with simulated and synthetic data, and then applied to fNIRS real data. The obtained results showed that the static version of our model is competitive against previous methodologies and demands a lower computational cost. Our model is more flexible than the previous methodologies by allowing the graphical structure to vary over time.

Functional concurrent, or varying-coefficient, regression models are a form of functional data analysis methods in which functional covariates and outcomes are collected concurrently. Two active areas of research for this class of models are identifying influential functional covariates and clustering their relations across observations. In various applications, researchers have applied and developed methods to address these objectives separately. However, no approach currently performs both tasks simultaneously. In this paper, we propose a fully Bayesian functional concurrent regression mixture model that simultaneously performs functional variable selection and clustering for subject-specific trajectories. Our approach introduces a novel spiked Ewens-Pitman attraction prior that identifies and clusters subjects’ trajectories marginally for each functional covariate while using similarities in subjects’ auxiliary covariate patterns to inform clustering allocation. Using simulated data, we evaluate the clustering, variable selection, and parameter estimation performance of our approach and compare its performance with alternative spiked processes. We then apply our method to functional data collected in a novel, smartphone-based smoking cessation intervention study to investigate individual-level dynamic relations between smoking behaviors and potential risk factors.

We present a flexible Bayesian semiparametric mixed model for longitudinal data analysis in the presence of potentially high-dimensional categorical covariates. Building on a novel hidden Markov tensor decomposition technique, our proposed method allows the fixed effects components to vary between dependent random partitions of the covariate space at different time points. The mechanism not only allows different sets of covariates to be included in the model at different time points but also allows the selected predictors’ influences to vary flexibly over time. Smooth time-varying additive random effects are used to capture subject specific heterogeneity. We establish posterior convergence guarantees for both function estimation and variable selection. We design a Markov chain Monte Carlo algorithm for posterior computation. We evaluate the method’s empirical performances through synthetic experiments and demonstrate its practical utility through real world applications.

Specification of the prior distribution for a Bayesian model is a central part of the Bayesian workflow for data analysis, but it is often difficult even for statistical experts. In principle, prior elicitation transforms domain knowledge of various kinds into well-defined prior distributions, and offers a solution to the prior specification problem. In practice, however, we are still fairly far from having usable prior elicitation tools that could significantly influence the way we build probabilistic models in academia and industry. We lack elicitation methods that integrate well into the Bayesian workflow and perform elicitation efficiently in terms of costs of time and effort. We even lack a comprehensive theoretical framework for understanding different facets of the prior elicitation problem.

Why are we not widely using prior elicitation? We analyse the state of the art by identifying a range of key aspects of prior knowledge elicitation, from properties of the modelling task and the nature of the priors to the form of interaction with the expert. The existing prior elicitation literature is reviewed and categorized in these terms. This allows recognizing under-studied directions in prior elicitation research, finally leading to a proposal of several new avenues to improve prior elicitation methodology.

Power-expected-posterior (PEP) methodology, which borrows ideas from the literature on power priors, expected-posterior priors and unit information priors, provides a systematic way to construct objective priors. The basic idea is to use imaginary training samples to update a (possibly improper) prior into a proper but minimally-informative one. In this work, we develop a novel definition of PEP priors for logistic regression models that relies on a Laplace expansion of the likelihood of the imaginary training sample. This approach has various advantages over previous proposals for non-informative priors in logistic regression, and can be easily extended to other generalized linear models. We study theoretical properties of the prior and provide a number of empirical studies that demonstrate superior performance both in terms of model selection and of parameter estimation, especially for heavy-tailed versions.

Completely random measures (CRMs) and their normalizations (NCRMs) offer flexible models in Bayesian nonparametrics. But their infinite dimensionality presents challenges for inference. Two popular finite approximations are truncated finite approximations (TFAs) and independent finite approximations (IFAs). While the former have been well-studied, IFAs lack similarly general bounds on approximation error, and there has been no systematic comparison between the two options. In the present work, we propose a general recipe to construct practical finite-dimensional approximations for homogeneous CRMs and NCRMs, in the presence or absence of power laws. We call our construction the automated independent finite approximation (AIFA). Relative to TFAs, we show that AIFAs facilitate more straightforward derivations and use of parallel computing in approximate inference. We upper bound the approximation error of AIFAs for a wide class of common CRMs and NCRMs — and thereby develop guidelines for choosing the approximation level. Our lower bounds in key cases suggest that our upper bounds are tight. We prove that, for worst-case choices of observation likelihoods, TFAs are more efficient than AIFAs. Conversely, we find that in real-data experiments with standard likelihoods, AIFAs and TFAs perform similarly. Moreover, we demonstrate that AIFAs can be used for hyperparameter estimation even when other potential IFA options struggle or do not apply.

In this paper, we study Bayesian asymptotic properties of the proportional hazards model where the link function is modeled by the generalized additive model. As the standard generalized additive model is, the generalized additive proportional hazards model is a useful tool in finding the nonlinearity of covariate effects to survival times. We develop a data-dependent sieve prior for the generalized additive link function and use the Bayesian bootstrap prior for the baseline hazard function. We prove that the posterior contraction rate of the generalized additive link function is minimax optimal up to a $logn$ term when the prior is carefully selected. By analyzing simulated as well as real data, we verify our theoretical results and compare with exisiting algorithms for the generalized additive proportional hazards model to illustrate that the proposed Bayesian model is a useful inference tool.

Bacteria reproduce clonally but most species recombine frequently, so that the ancestral process is best captured using an ancestral recombination graph. This graph model is often too complex to be used in an inferential setup, but it can be approximated for example by the ClonalOrigin model. Inference in the ClonalOrigin model is performed via a Reversible-Jump Markov Chain Monte Carlo algorithm, which attempts to jointly explore: the recombination rate, the number of recombination events, the departure and arrival points on the clonal genealogy for each recombination event, and the range of genomic sites affected by each recombination event. However, the Reversible-Jump algorithm often performs poorly due to the complexity of the target distribution since it needs to explore spaces of different dimensions. Recent developments in Bayesian computation methodology have provided ways to improve existing methods and code, but are not well-known outside the statistics community. We show how exploiting one of these new computational methods can lead to faster inference under the ClonalOrigin model.

In this paper, we propose a new Bayesian inference method for a high-dimensional sparse factor model that allows both the factor dimensionality and the sparse structure of the loading matrix to be inferred. The novelty is to introduce a certain dependence between the sparsity level and the factor dimensionality, which leads to adaptive posterior concentration while keeping computational tractability. We show that the posterior distribution asymptotically concentrates on the true factor dimensionality, and more importantly, this posterior consistency is adaptive to the sparsity level of the true loading matrix and the noise variance. We also prove that the proposed Bayesian model attains the optimal detection rate of the factor dimensionality in a more general situation than those found in the literature. Moreover, we obtain a near-optimal posterior concentration rate of the covariance matrix. Numerical studies are conducted and show the superiority of the proposed method compared with other competitors.

The prominent Bernstein – von Mises (BvM) Theorem claims a kind of approximation of the posterior distribution by a Gaussian one with the covariance close to the inverse of the total Fisher information matrix. A more general Laplace approximation result states a similar Gaussian approximation of the posterior with the parameters depending on the prior. These two results build a basis for Bayesian inference and uncertainty quantification in a rather general situation. Spokoiny and Panov (2021) offered a new look at this problem which allows to state rather strong results on the quality of Gaussian approximation in non-asymptotic and dimension free form assuming linearity and concavity of log-likelihood function which can be misspecified. The established results provide explicit non-asymptotic bounds on the quality of a Gaussian approximation of the posterior distribution in total variation distance in terms of the so called effective dimension ${\mathtt{p}}_G$ defined as interplay between information contained in the data and in the prior distribution. This paper substantially improves and further develops the results from Spokoiny and Panov (2021) using the recent progress on high dimensional Laplace approximation. We address the question of effective and critical dimension in Bayesian inference, the relations between Laplace approximation and Bernstein–von Mises Theorem, and, particularly, the use of posterior mean instead of Maximum A Posteriori Probability estimator in Bayesian inference. The results are illustrated for the case of log-density estimation.

Dynamic Linear Models (DLMs) are commonly employed for time series analysis due to their versatile structure, simple recursive updating, ability to handle missing data, and probabilistic forecasting. However, the options for count time series are limited: Gaussian DLMs require continuous data, while Poisson-based alternatives often lack sufficient modeling flexibility. We introduce a novel semiparametric methodology for count time series by warping a Gaussian DLM. The warping function has two components: a (nonparametric) transformation operator that provides distributional flexibility and a rounding operator that ensures the correct support for the discrete data-generating process. We develop conjugate inference for the warped DLM, which enables analytic and recursive updates for the state space filtering and smoothing distributions. We leverage these results to produce customized and efficient algorithms for inference and forecasting, including Monte Carlo simulation for offline analysis and an optimal particle filter for online inference. This framework unifies and extends a variety of discrete time series models and is valid for natural counts, rounded values, and multivariate observations. Simulation studies illustrate the excellent forecasting capabilities of the warped DLM. The proposed approach is applied to a multivariate time series of daily overdose counts and demonstrates both modeling and computational successes.

Copula-based models provide a great deal of flexibility in modelling multivariate distributions, allowing for the specifications of models for the marginal distributions separately from the dependence structure (copula) that links them to form a joint distribution. Choosing a class of copula models is not a trivial task and its misspecification can lead to wrong conclusions. We introduce a novel class of grid-uniform copula functions, which is dense in the space of all continuous copula functions in a Hellinger sense. We propose a Bayesian model based on this class which posterior distribution is strongly consistent and develop an automatic Markov chain Monte Carlo algorithm for exploring the corresponding posterior distribution. The methodology is illustrated by means of simulated data and compared to the main existing approach.

State space models (SSMs) provide a flexible framework for modeling complex time series via a latent stochastic process. Inference for nonlinear, non-Gaussian SSMs is often tackled with particle methods that do not scale well to long time series. The challenge is two-fold: not only do computations scale linearly with time, as in the linear case, but particle filters additionally suffer from increasing particle degeneracy with longer series. Stochastic gradient MCMC methods have been developed to scale Bayesian inference for finite-state hidden Markov models and linear SSMs using buffered stochastic gradient estimates to account for temporal dependencies. We extend these stochastic gradient estimators to nonlinear SSMs using particle methods. We present error bounds that account for both buffering error and particle error in the case of nonlinear SSMs that are log-concave in the latent process. We evaluate our proposed particle buffered stochastic gradient using stochastic gradient MCMC for inference on both long sequential synthetic and minute-resolution financial returns data, demonstrating the importance of this class of methods.

In many common situations, a Bayesian credible interval will be, given the same data, very similar to a frequentist confidence interval, and researchers will interpret these intervals in a similar fashion. However, no predictable similarity exists when credible intervals are based on model-averaged posteriors whenever one of the two nested models under consideration is a so called “point-null”. Not only can this model-averaged credible interval be quite different than the frequentist confidence interval, in some cases it may be undefined. This is a lesser-known correlate of the Jeffreys-Lindley paradox and is of particular interest given the popularity of the Bayes factor for testing point-null hypotheses.

Deterministic compartmental models are predominantly used in the modeling of infectious diseases, though stochastic models are considered more realistic, yet are complicated to estimate due to missing data. In this paper we present a novel algorithm for estimating the stochastic Susceptible-Infected-Recovered and Susceptible-Exposed-Infected-Recovered (SIR/SEIR) epidemic model within a Bayesian framework, which can be readily extended to more complex stochastic compartmental models. Specifically, based on the infinitesimal conditional independence properties of the model, we are able to find a proposal distribution for a Metropolis algorithm which is very close to the correct posterior distribution. In particular, it acts as a very good proposal for the unknown number of events, such as the number of infected individuals. Finding such a good proposal for this has historically been a difficult problem. As a consequence, rather than perform a Metropolis step updating one missing data point at a time, we are able to have a single proposal for the entire set of missing observations. A number of real data illustrations and the necessary mathematical theory supporting our results are presented.

Diagnosing convergence of Markov chain Monte Carlo is crucial and remains an essentially unsolved problem. Among the most popular methods, the potential scale reduction factor, commonly named $\hat{R}$, is an indicator that monitors the convergence of output chains to a target distribution, based on a comparison of the between- and within-variances. Several improvements have been suggested since its introduction in the 90s. Here, we aim at better understanding the $\hat{R}$ behavior by proposing a localized version that focuses on quantiles of the target distribution. This new version relies on key theoretical properties of the associated population value. It naturally leads to proposing a new indicator ${\hat{R}}_{\infty }$, which is shown to allow both for localizing the Markov chain Monte Carlo convergence in different quantiles of the target distribution, and at the same time for handling some convergence issues not detected by other $\hat{R}$ versions.

Standard approaches for functional principal components analysis rely on an eigendecomposition of a smoothed covariance surface in order to extract the orthonormal eigenfunctions representing the major modes of variation in a set of functional data. This approach can be a computationally intensive procedure, especially in the presence of large datasets with irregular observations. In this article, we develop a variational Bayesian approach, which aims to determine the Karhunen-Loève decomposition directly without smoothing and estimating a covariance surface. More specifically, we incorporate the notion of variational message passing over a factor graph because it removes the need for rederiving approximate posterior density functions if there is a change in the model. Instead, model changes are handled by changing specific computational units, known as fragments, within the factor graph – we demonstrate this with an extension to multilevel functional data. Indeed, this is the first article to address a functional data model via variational message passing. Our approach introduces three new fragments that are necessary for Bayesian functional principal components analysis. We present the computational details, a set of simulations for assessing the accuracy and speed of the variational message passing algorithm and an application to United States temperature data.

Quantile feature selection over correlated multivariate time series data has always been a methodological challenge and is an open problem. In this paper, we propose a general Bayesian dimension reduction methodology for feature selection in high-dimensional joint quantile time series analysis, under the name of the quantile feature selection time series (QFSTS) model. The QFSTS model is a general structural time series model, where each component yields an additive contribution to the time series modeling with direct interpretations. Its flexibility is compound in the sense that users can add/deduct components for each time series and each time series can have its own specific valued components of different sizes. Feature selection is conducted in the quantile regression component, where each time series has its own pool of contemporaneous external predictors allowing nowcasting. Bayesian methodology in extending feature selection to the quantile time series research area is developed using multivariate asymmetric Laplace distribution, spike-and-slab prior setup, the Metropolis-Hastings algorithm, and the Bayesian model averaging technique, all implemented consistently in the Bayesian paradigm. The QFSTS model requires small datasets to train and converges fast. Extensive examinations confirmed that the QFSTS model has superior performance in feature selection, parameter estimation, and forecast.

Estimating heterogeneous treatment effects has drawn increasing attention in medical studies, considering that patients with divergent features can undergo a different progression of disease even with identical treatment. Such heterogeneity can co-occur with a cured fraction for biomedical studies with a time-to-event outcome and further complicates the quantification of treatment effects. This study considers a joint framework of Bayesian causal forest and accelerated failure time cure model to capture the cured proportion and treatment effect heterogeneity through three separate Bayesian additive regression trees. Under the potential outcomes framework, conditional and sample average treatment effects within the uncured subgroup are derived on the scale of log survival time subject to right-censoring, and treatment effects on the scale of survival probability are derived for each individual. Bayesian backfitting Markov chain Monte Carlo algorithm with the Gibbs sampler is conducted to estimate the causal effects. Simulation studies show the satisfactory performance of the proposed method. The proposed model is then applied to a breast cancer dataset extracted from the SEER database to demonstrate its usage in detecting heterogeneous treatment effects and cured subgroups. Combined with popular mitigation strategies, the proposed method can also alleviate confounding induced by immortal time bias.

Differential equations based on physical principals are used to represent complex dynamic systems in all fields of science and engineering. When known, these equations have been shown to well represent real-world dynamics. However, since the true dynamics of complex systems are generally unknown, learning the governing equations can improve our understanding of the mechanisms driving the systems. Here, we develop a Bayesian approach to data-driven discovery of nonlinear spatio-temporal dynamic equations. Our approach can accommodate measurement error and missing data, both of which are common in real-world data, and accounts for parameter uncertainty. The proposed framework is illustrated using three simulated systems with varying amounts of measurement uncertainty and missing data and applied to a real-world system to infer the temporal evolution of the vorticity of the streamfunction.

Interactions between actors are frequently represented using a network. The latent position model is widely used for analysing network data, whereby each actor is positioned in a latent space. Inferring the dimension of this space is challenging. Often, for simplicity, two dimensions are used or model selection criteria are employed to select the dimension, but this requires choosing a criterion and the computational expense of fitting multiple models. Here the latent shrinkage position model (LSPM) is proposed which intrinsically infers the effective dimension of the latent space. The LSPM employs a Bayesian nonparametric multiplicative truncated gamma process prior that ensures shrinkage of the variance of the latent positions across higher dimensions. Dimensions with non-negligible variance are deemed most useful to describe the observed network, inducing automatic inference on the latent space dimension. While the LSPM is applicable to many network types, logistic and Poisson LSPMs are developed here for binary and count networks respectively. Inference proceeds via a Markov chain Monte Carlo algorithm, where novel surrogate proposal distributions reduce the computational burden. The LSPM’s properties are assessed through simulation studies, and its utility is illustrated through application to real network datasets. Open source software assists wider implementation of the LSPM.

Markov chain Monte Carlo methods for exponential family models with intractable normalizing constant, such as the exchange algorithm, require simulations of the sufficient statistics at every iteration of the Markov chain, which often result in expensive computations. Surrogate models for the likelihood function have been developed to accelerate inference algorithms in this context. However, these surrogate models tend to be relatively inflexible, and often provide a poor approximation to the true likelihood function. In this article, we propose the use of a warped, gradient-enhanced, Gaussian process surrogate model for the likelihood function, which jointly models the sample means and variances of the sufficient statistics, and uses warping functions to capture covariance nonstationarity in the input parameter space. We show that both the consideration of nonstationarity and the inclusion of gradient information can be leveraged to obtain a surrogate model that outperforms the conventional stationary Gaussian process surrogate model when making inference, particularly in regions where the likelihood function exhibits a phase transition. We also show that the proposed surrogate model can be used to improve the effective sample size per unit time when embedded in exact inferential algorithms. The utility of our approach in speeding up inferential algorithms is demonstrated on simulated and real-world data.

Simulation-based calibration checking (SBC) is a practical method to validate computationally-derived posterior distributions or their approximations. In this paper, we introduce a new variant of SBC to alleviate several known problems. Our variant allows the user to in principle detect any possible issue with the posterior, while previously reported implementations could never detect large classes of problems including when the posterior is equal to the prior. This is made possible by including additional data-dependent test quantities when running SBC. We argue and demonstrate that the joint likelihood of the data is an especially useful test quantity. Some other types of test quantities and their theoretical and practical benefits are also investigated. We provide theoretical analysis of SBC, thereby providing a more complete understanding of the underlying statistical mechanisms. We also bring attention to a relatively common mistake in the literature and clarify the difference between SBC and checks based on the data-averaged posterior. We support our recommendations with numerical case studies on a multivariate normal example and a case study in implementing an ordered simplex data type for use with Hamiltonian Monte Carlo. The SBC variant introduced in this paper is implemented in the SBC R package.

DNA methylation datasets in cancer studies are comprised of measurements on a large number of genomic locations called cytosine-phosphate-guanine (CpG) sites with complex correlation structures. A fundamental goal of these studies is the development of statistical techniques that can identify disease genomic signatures across multiple patient groups defined by different experimental or biological conditions. We propose BayesDiff, a nonparametric Bayesian approach for differential analysis relying on a novel class of first order mixture models called the Sticky Pitman-Yor process or two-restaurant two-cuisine franchise (2R2CF). The BayesDiff methodology flexibly utilizes information from all CpG sites or biomarker probes, adaptively accommodates any serial dependence due to the widely varying inter-probe distances, and makes posterior inferences about the differential genomic signature of patient groups. Using simulation studies, we demonstrate the effectiveness of the BayesDiff procedure relative to existing statistical techniques for differential DNA methylation. The methodology is applied to analyze a gastrointestinal (GI) cancer dataset exhibiting serial correlation and complex interaction patterns. The results support and complement known aspects of DNA methylation and gene association in upper GI cancers.

We consider the problem of performing inference for a continuous treatment effect on a binary outcome variable while controlling for high dimensional baseline covariates. We propose a novel Bayesian framework for performing inference for the desired low-dimensional parameter in a high-dimensional logistic regression model. While it is relatively easier to address this problem in linear regression, the nonlinearity of the logistic regression poses additional challenges that make it difficult to orthogonalize the effect of the treatment variable from the nuisance variables. Our proposed approach provides the first Bayesian alternative to the recent frequentist developments and can incorporate available prior information on the parameters of interest, which plays a crucial role in practical applications. In addition, the proposed approach incorporates uncertainty in orthogonalization in high dimensions instead of relying on a single instance of orthogonalization as done by frequentist methods. We provide uniform convergence results that show the validity of credible intervals resulting from the posterior. Our method has competitive empirical performance when compared with state-of-the-art methods.

Normalized random measures with independent increments (NRMIs) represent a large class of Bayesian nonparametric priors and are widely used in the Bayesian nonparametric framework. In this paper, we provide the posterior consistency analysis for these NRMIs through their characterizing Lévy intensities. Assumptions are introduced on the Lévy intensities to analyse the posterior consistency and are verified with multiple interesting examples. Another focus of the paper is the Bernstein-von Mises theorem for a particular subclass of NRMIs, namely the normalized generalized gamma processes (NGGP). When the Bernstein-von Mises theorem is applied to construct credible sets, in addition to the usual form, there will be an additional bias term on the left endpoint closely related to the number of atoms of the true distribution in the discrete case. We also discuss the effect of the estimators for the model parameters of the NGGP under the Bernstein-von Mises convergence. Finally, to further illustrate the impact of the bias correction term in the construction of credible sets, we present a numerical example to demonstrate numerically how the bias correction affects the coverage of the true value.

Bayesian cross-validation (CV) is a popular method for predictive model assessment that is simple to implement and broadly applicable. A wide range of CV schemes is available for time series applications, including generic leave-one-out (LOO) and K-fold methods, as well as specialized approaches intended to deal with serial dependence such as leave-future-out (LFO), h-block, and $hv$-block.

Existing large-sample results show that both specialized and generic methods are applicable to models of serially-dependent data. However, large sample consistency results overlook the impact of sampling variability on accuracy in finite samples. Moreover, the accuracy of a CV scheme depends on many aspects of the procedure. We show that poor design choices can lead to elevated rates of adverse selection.

In this paper, we consider the problem of identifying the regression component of an important class of models of data with serial dependence, autoregressions of order p with q exogenous regressors ($\mathsf{ARX}(p,q)$), under the logarithmic scoring rule. We show that when serial dependence is present, scores computed using the joint (multivariate) density have lower variance and better model selection accuracy than the popular pointwise estimator. In addition, we present a detailed case study of the special case of $\mathsf{ARX}$ models with fixed autoregressive structure and variance. For this class, we derive the finite-sample distribution of the CV estimators and the model selection statistic. We conclude with recommendations for practitioners.

Mobile health (mHealth) methods allow researchers to monitor study participants in their natural environments in order to improve health-related outcomes through behavior change. MHealth investigators are interested in understanding the feasibility of supplementing or even replacing actively collected data with passively collected data to reduce participant burden, while also increasing the temporal resolution of the data. In this work, we propose a novel Bayesian dynamic functional variable selection method to explore the relations between multimodal mHealth data collected on different time scales. Specifically, our approach leverages spiked hierarchical species sampling priors to identify critical moments when a participant experiences momentary spikes in a low-frequency outcome, characterize high-frequency outcome trajectories which are related to these critical moments, and cluster the relational trends to explore potential subpopulations. We introduce continuous-time multistate Markov model priors to inform selection based on information learned at previous assessments. We demonstrate the variable selection and clustering performance of our model in various simulation settings motivated by the data structures found in mHealth studies. We then apply our model to multimodal intensive longitudinal data collected in the Pathways between Socioeconomic Status and Behavioral Cancer Risk Factors Study to explore relations between physical activity passively collected with accelerometers and mood actively collected with ecological momentary assessment methods.

In environmental health research there is often interest in the effect of an exposure on a health outcome assessed on the same day and several subsequent days or lags. Distributed lag nonlinear models (DLNM) are a well-established statistical framework for estimating an exposure-lag-response function. We propose methods to allow for prior information to be incorporated into DLNMs. First, we impose a monotonicity constraint in the exposure-response at lagged time periods which matches with knowledge on how biological mechanisms respond to increased levels of exposures. Second, we introduce variable selection into the DLNM to identify lagged periods of susceptibility with respect to the outcome of interest. The variable selection approach allows for direct application of informative priors on which lags have nonzero association with the outcome. We propose a tree-of-trees model that uses two layers of trees: one for splitting the exposure time frame and one for fitting exposure-response functions over different time periods. We introduce a zero-inflated alternative to the tree splitting prior in Bayesian additive regression trees to allow for lag selection and the addition of informative priors. We develop a computational approach for efficient posterior sampling and perform a comprehensive simulation study to compare our method to existing DLNM approaches. We apply our method to estimate time-lagged extreme temperature relationships with mortality during summer or winter in Chicago, IL.

It is common to hold prior beliefs that are not characterized by points in the parameter space but instead are relational in nature and can be described by a linear subspace. While some previous work has been done to account for such prior beliefs, the focus has primarily been on point estimators within a regression framework. We argue, however, that prior beliefs about parameters ought to be encoded into the prior distribution rather than in the formation of a point estimator. In this way, the prior beliefs help shape all inference. Through exponential tilting, we propose a fully generalizable method of taking existing prior information from, e.g., a pilot study, and combining it with additional prior beliefs represented by parameters lying on a linear subspace. We provide computationally efficient algorithms for posterior inference that, once inference is made using a non-tilted prior, does not depend on the sample size. We illustrate our proposed approach on an antihypertensive clinical trial dataset where we shrink towards a power law dose-response relationship, and on monthly influenza and pneumonia data where we shrink moving average lag parameters towards smoothness. Software to implement the proposed approach is provided in the R package SUBSET available on GitHub.

Bayesian nonparametric mixture models offer a rich framework for model-based clustering. We consider the situation where the kernel of the mixture is available only up to an intractable normalizing constant. In this case, the most commonly used Markov chain Monte Carlo (MCMC) methods are unsuitable. We propose an approximate Bayesian computational (ABC) strategy, whereby we approximate the posterior to avoid the intractability of the kernel. We derive an ABC-MCMC algorithm which combines (i) the use of the predictive distribution induced by the nonparametric prior as proposal and (ii) the use of the Wasserstein distance and its connection to optimal matching problems. To overcome the sensibility concerning the parameters of our algorithm, we further propose an adaptive strategy. We illustrate the use of the proposed algorithm with several simulation studies and an application on real data, where we cluster a population of networks, comparing its performance with standard MCMC algorithms and validating the adaptive strategy.

Species-sampling problems (SSPs) refer to a vast class of statistical problems calling for the estimation of (discrete) functionals of the unknown species composition of an unobservable population. A common feature of SSPs is their invariance with respect to species labeling, which is at the core of the Bayesian nonparametric (BNP) approach to SSPs under the popular Pitman-Yor process (PYP) prior. In this paper, we develop a BNP approach to SSPs that are not “invariant” to species labeling, in the sense that an ordering or ranking is assigned to species’ labels. Inspired by the population genetics literature on age-ordered alleles’ compositions, we study the following SSP with ordering: given an observable sample from an unknown population of individuals belonging to species (alleles), with species’ labels being ordered according to weights (ages), estimate the frequencies of the first r order species’ labels in an enlarged sample obtained by including additional unobservable samples. By relying on an ordered PYP prior, we obtain an explicit posterior distribution of the first r order frequencies, with estimates being of easy implementation and computationally efficient. We apply our approach to the analysis of genetic variation, showing its effectiveness in estimating the frequency of the oldest allele, and then we discuss other potential applications.

The increasing prevalence of network data in a vast variety of fields and the need to extract useful information out of them have spurred fast developments in related models and algorithms. Among the various learning tasks with network data, community detection, the discovery of node clusters or “communities,” has arguably received the most attention in the scientific community. In many real-world applications, the network data often come with additional information in the form of node or edge covariates that should ideally be leveraged for inference. In this paper, we add to a limited literature on community detection for networks with covariates by proposing a Bayesian stochastic block model with a covariate-dependent random partition prior. Under our prior, the covariates are explicitly expressed in specifying the prior distribution on the cluster membership. Our model has the flexibility of modeling uncertainties of all the parameter estimates including the community membership. Importantly, and unlike the majority of existing methods, our model has the ability to learn the number of the communities via posterior inference without having to assume it to be known. Our model can be applied to community detection in both dense and sparse networks, with both categorical and continuous covariates, and our MCMC algorithm is very efficient with good mixing properties. We demonstrate the superior performance of our model over existing models in a comprehensive simulation study and an application to two real datasets.

Dirichlet processes and their extensions have reached a great popularity in Bayesian nonparametric statistics. They have also been introduced for spatial and spatio-temporal data, as a tool to analyze and predict surfaces. A popular approach to Dirichlet processes in a spatial setting relies on a stick-breaking representation of the process, where the dependence over space is described in the definition of the stick-breaking probabilities. Extensions to include temporal dependence are still limited, however it is important, in particular for those phenomena which may change rapidly over time and space, with many local changes. In this work, we propose a Dirichlet process where the stick-breaking probabilities are defined to incorporate both spatial and temporal dependence. We will show that this approach is not a simple extension of available methodologies and can outperform available approaches in terms of prediction accuracy. The approach results in a predictive model that does not rely on the assumptions of Gaussianity or the separability of time and space.

In this paper, we study the finite approximation of the completely random measure (CRM) by truncating its Ferguson-Klass representation. The approximation is obtained by keeping the N largest atom weights of the CRM unchanged and combining the smaller atom weights into a single term. We develop the simulation algorithms for the approximation and characterise its posterior distribution, for which a blocked Gibbs sampler is devised. We demonstrate the usage of the approximation in two models. The first assumes such an approximation as the mixing distribution of a Bayesian nonparametric mixture model and leads to a finite approximation to the model posterior. The second concerns the finite approximation to the Caron-Fox model. Examples and numerical implementations are given based on the gamma, stable and generalised gamma processes.

Heavy tails are often found in practice, and yet they are an Achilles heel of a variety of mainstream random probability measures such as the Dirichlet process (DP). The first contribution of this paper focuses on characterizing the tails of the so-called normalized generalized gamma (NGG) process. We show that the right tail of an NGG process is heavy-tailed provided that the centering distribution is itself heavy-tailed; the DP is the only member of the NGG class that fails to obey this convenient property. A second contribution of the paper rests on the development of two classes of heavy-tailed mixture models and the assessment of their relative merits. Multivariate extensions of the proposed heavy-tailed mixtures are devised here, along with a predictor-dependent version, to learn about the effect of covariates on a multivariate heavy-tailed response. The simulation study suggests that the proposed method performs well in various scenarios, and we showcase the application of the proposed methods in a neuroscience dataset.

In recent years, new cancer treatments have improved survival in multiple histologies. Some of these therapeutics, and in particular treatment combinations, are often associated with severe treatment-related adverse events (AEs). Therefore, It is important to identify alternative de-intensified therapies, such as dose-reduced therapies with reduced AEs and similar efficacy. We introduce a sequential design for multi-arm de-intensification studies. Based on joint modeling of toxicity and efficacy endpoints, the design evaluates multiple de-intensified therapies at different dose levels, one at a time. We study the utility of the design in oropharynx cancer de-intensification studies. We use a joint Bayesian model for efficacy and toxicity outcomes to define decision rules at interim and final analyses. Interim decisions include early termination of the study due to inferior survival of experimental arms compared to a standard of care (SOC), and the transitions from one de-intensified treatment arm to another with a further reduced dose when there is sufficient evidence of non-inferior survival. We evaluate the operating characteristics of the design using data from recent de-intensification studies in oropharynx cancer.

There has been increased research interest in the subfield of sparse Bayesian factor analysis with shrinkage priors, which achieve additional sparsity beyond the natural parsimonity of factor models. In this spirit, we estimate the number of common factors in the widely applied sparse latent factor model with spike-and-slab priors on the factor loadings matrix. Our framework leads to a natural, efficient and simultaneous coupling of model estimation and selection on one hand and model identification and rank estimation (number of factors) on the other hand. More precisely, by embedding the unordered generalised lower triangular loadings representation into overfitting sparse factor modelling, we obtain posterior summaries regarding factor loadings, common factors as well as the factor dimension via postprocessing draws from our efficient and customized Markov chain Monte Carlo scheme.

Clinical trials often allow patients in the control arm to switch to the treatment arm if their physical conditions are worse than certain tolerance levels. For instance, treatment switching arises in the Concorde clinical trial, which aims to assess causal effects on the time-to-disease progression or death of immediate versus deferred treatment with zidovudine among patients with asymptomatic HIV infection. The Intention-To-Treat analysis does not measure the effect of the actual receipt of the treatment and ignores the information on treatment switching. Other existing methods reconstruct the outcome a patient would have had if they had not switched under strong assumptions. Departing from the literature, we re-define the problem of treatment switching using principal stratification and focus on causal effects for patients belonging to subpopulations defined by the switching behavior under control. We use a Bayesian approach to inference, taking into account that $(i)$ switching happens in continuous time; $(ii)$ switching time is not defined for patients who never switch in a particular experiment; and $(iii)$ survival time and switching time are subject to censoring. We apply this framework to analyze the synthetic data based on the Concorde study.

Additive spatial statistical models with weakly stationary process assumptions have become standard in spatial statistics. However, one disadvantage of such models is the computation time, which rapidly increases with the number of data points. The goal of this article is to apply an existing subsampling strategy to standard spatial additive models and to derive the spatial statistical properties. We call this strategy the “spatial data subset model” (SDSM) approach, which can be applied to big datasets in a computationally feasible way. Our approach has the advantage that one does not require any additional restrictive model assumptions. That is, computational gains increase as model assumptions are removed when using our model framework. This provides one solution to the computational bottlenecks that occur when applying methods such as Kriging to “big data”. We provide several properties of this new spatial data subset model approach in terms of moments, sill, nugget, and range under several sampling designs. An advantage of our approach is that it subsamples without throwing away data, and can be implemented using datasets of any size that can be stored. We present the results of the spatial data subset model approach on simulated datasets, and on a large dataset consists of 150,000 observations of daytime land surface temperatures measured by the MODIS instrument onboard the Terra satellite.

As a regression technique in spatial statistics, the spatiotemporally varying coefficient model (STVC) is an important tool for discovering nonstationary and interpretable response-covariate associations over both space and time. However, it is difficult to apply STVC for large-scale spatiotemporal analyses due to its high computational cost. To address this challenge, we summarize the spatiotemporally varying coefficients using a third-order tensor structure and propose to reformulate the spatiotemporally varying coefficient model as a special low-rank tensor regression problem. The low-rank decomposition can effectively model the global patterns of large data sets with a substantially reduced number of parameters. To further incorporate the local spatiotemporal dependencies, we use Gaussian process (GP) priors on the spatial and temporal factor matrices. We refer to the overall framework as Bayesian Kernelized Tensor Regression (BKTR), and kernelized tensor factorization can be considered a new and scalable approach to modeling multivariate spatiotemporal processes with a low-rank covariance structure. For model inference, we develop an efficient Markov chain Monte Carlo (MCMC) algorithm, which uses Gibbs sampling to update factor matrices and slice sampling to update kernel hyperparameters. We conduct extensive experiments on both synthetic and real-world data sets, and our results confirm the superior performance and efficiency of BKTR for model estimation and parameter inference.

Within the field of record linkage, Bayesian methods have the crucial advantage of quantifying uncertainty from imperfect linkages. However, current implementations of Bayesian Fellegi-Sunter models are computationally intensive, making them challenging to use on larger-scale record linkage tasks. To address these computational difficulties, we propose fast beta linkage (fabl), an extension to the Beta Record Linkage (BRL) method of Sadinle (2017). Specifically, we use independent prior distributions over the matching space, allowing us to use hashing techniques that reduce computational overhead. This also allows us to complete pairwise record comparisons over large data files through parallel computing and to reduce memory costs through a new technique called storage efficient indexing. Through simulations and two case studies, we show that fabl can have markedly increased speed with minimal loss of accuracy when compared to BRL.

Due to their flexibility Gaussian processes are a well-known Bayesian framework for nonparametric function estimation. Integrating inequality constraints, such as monotonicity, convexity, and boundedness, into Gaussian process models significantly improves prediction accuracy and yields more realistic credible intervals in various real-world data applications. The finite-dimensional Gaussian process approximation, originally proposed in Maatouk and Bay (2017) is considered. This method involves approximating a parent GP by utilizing a finite-dimensional GP obtained through appropriate basis expansions. It satisfies interpolation conditions and handles a wide range of inequality constraints everywhere. Our contribution in this paper is threefold. First, we extend this approach to handle noisy observations and multiple, more general convex and non-convex constraints. Second, we propose new basis functions in order to extend the smoothness of sample paths to differentiability of class $C^p$, for any $p\ge 1$. Third, we examine its behavior in specific scenarios such as monotonicity with flat regions and boundedness near lower and/or upper bounds. In that case, we show that, unlike the Maximum a posteriori (MAP) estimate, the mean a posteriori (mAP) estimate fails to capture flat regions. To address this issue, we propose incorporating multiple constraints, such as monotonicity with bounded slope constraints. According to the theoretical convergence and based on a variety of numerical experiments, the MAP estimate behaves well and outperforms the mAP estimate in terms of prediction accuracy. The performance of the proposed approach is confirmed through synthetic and real-world data studies.