Introductory Chapter: Development of Bayesian Inference

1. Introduction

Bayesian inference derives from Bayesian theory, which depicts the probability of occurrence of an event given some prior information. Due to the huge advances in computational and modeling techniques, Bayesian inference has been increasingly become an important tool for data analysis in the Bayesian framework and has widely been applied to various fields, including social science, engineering, philosophy, medicine, sport, law and psychology, for parameter/nonparameter estimation, hypothesis test, and prediction. Various Bayesian methods including Markov chain Monte Carlo, objective Bayesian method, subjective Bayesian method, approximate Bayesian computation, and variational Bayesian methods have been developed to make Bayesian inference on various problems such as large-scale image classification and cluster analysis of microarray, and models including parametric, nonparametric, semiparametric models, and other complicated models such as joint models of survival data and longitudinal data, graphical models, computer models, neural network models, and spatial econometric models. In particular, in the big data era, various Bayesian fronts including theories, methods, and computational algorithms have been developed for accommodating the applications of AI and data science in recent years [1], for example, the prior learning, Bayes factor evaluation, Bayesian variable selection, robust Bayesian inference, variational Bayesian inference, resampling, approximation of posterior distribution, approximate Bayesian computation, and debias methods for high−/ultrahigh-dimensional data, multisource heterogeneous data, imbalanced data, missing data, and data stream. But there are some challenging problems, for example, how to balance the computational times and statistical efficiency, design efficient Bayesian computational algorithm and robust sampling schemes for big/massive data, distributed data and streaming data in the privacy protection and the defense of malicious attacks framework, and modeling so that they can adapt to the development of AI and the requirement of data mining, to be addressed and solved for Bayesian inference. In what follows, we will introduce the recent development and some topics of interest on Bayesian inference.

2. Bayesian estimation

For statistical models, Bayesian estimation is usually obtained from its posterior distribution based on Bayes theory. In general, Bayesian estimation includes Bayesian estimation of parameters and nonparametric functions. For parametric Bayesian estimation, we need first to specify the prior distribution of the parameter and then evaluate its posterior mean/median (i.e., Bayesian estimation of parameter) from its posterior distributions if the quadratic/mean-absolute loss function is used. For nonparametric Bayesian estimation, we need first approximate nonparametric function via some proper method such as B-splines and P-splines [2], i.e., parameterized approximation to nonparametric function, which leads to a parameterized model, and then employ Bayesian idea of parametric models to evaluate Bayesian estimates of parameter and nonparametric function. In what follows, we introduce how to evaluate Bayesian estimates of parameters or nonparametric functions in a relatively complicated model (e.g., random effects model/latent variable model).

In a latent variable model with missing response data, we assume the following form:

$$Y_i=X_i\beta +\mathrm{\Lambda }{\omega }_i+{\epsilon }_i,{\eta }_i=\mathrm{\Pi }{\eta }_i+\mathrm{\Gamma }{\xi }_i+{ϵ}_i,i=1,2,\dots ,n,$$E1

where $Y_i$ is a $p\times 1$ vector of manifest variables including continuous and categorical variables [3], $X_i$ is a $p\times q$ matrix of covariates, ${\omega }_i$ is a $r\times 1$ vector of latent variables, ${ϵ}_i$ and ${\epsilon }_i$ are the $p\times 1$ and $r_1\times 1$ vectors of measurement errors, respectively, $\beta$ is a $q\times 1$ vector of unknown parameters, $\mathrm{\Lambda }$ is a $p\times r$ factor loading matrix, ${\omega }_i={\left({\eta }_i^{\mathit{T}},{\xi }_i^{\mathit{T}}\right)}^{\mathit{T}}$ in which ${\eta }_i$ and ${\xi }_i$ are the $r_1\times 1$ and $r_2\times 1$ sub-vectors of ${\omega }_i$ and $r_1+r_2=r$, respectively, $\mathrm{\Pi }$ and $\mathrm{\Gamma }$ are $r_1\times r_1$ and $r_1\times r_2$ matrices of unknown parameters. It is also assumed that $\mid I-\mathrm{\Pi }\mid \ne 0$, ${\xi }_i$ follows a normal distribution or an unknown distribution, and $y_{\mathit{ij}}$‘s are subject to missingness, where $I$ is a $r_1\times r_1$ identity matrix, $y_{\mathit{ij}}$ is the $j$th component of $Y_i$ for $i=1,\dots ,n,j=1,\dots ,p$.

In general, a simple and standard assumption for the distributions of ${\epsilon }_i$, ${ϵ}_i$_, and ${\xi }_i$ is to follow some parametric density family such as skew-normal/skew-t/skew-normal-cauchy/skew-symmetric-Laplace distribution [2] or exponential family distribution [3] or normal distribution [4], or unequal time autoregression series or their mixture. But this assumption may be unreasonable or too restrictive. To this end, some alternative methods have been developed to relax their parametric distribution assumption. For example, let ${\epsilon }_i$ or ${\xi }_i$ follow an unknown distribution, which is approximated by a Dirichlet process prior [2], spiked Dirichlet process prior [5], truncated centered Dirichlet process prior [6], or a class of smooth densities that are approximated by the semiparametric approach of Gallant and Nychka [7], or unknown distribution with its quantiles specified leading to the well-known quantile regression models, or a Bayesian neural network approach to learning unknown distribution [8].

To introduce missing data, let ${\delta }_{\mathit{ij}}$ be an indicator of missing value $y_{\mathit{ij}}$, i.e., ${\delta }_{\mathit{ij}}=1$ if $y_{\mathit{ij}}$ is observed, and ${\delta }_{\mathit{ij}}=0$ if $y_{\mathit{ij}}$ is missing. In this case, we usually need to assume a missingness data mechanism, for example, missing at completely random or missing at random or missing not at random (also called nonignorable missing), and then specify a parametric or semiparmetric model for the considered missingness data mechanism. For example, see Lee and Tang [4] for the logistic regression model, Kim and Yu [9] and Tang, Zhao, and Zhu [10] for the exponential tilting model, and Wang and Tang [11] for the probit regression model together with latent variables. Also, one can train a missingness data mechanism model via a data-driven machine learning method [12], which is a completely new and not yet studied topic, including its implementation, algorithm, and theories. We are working on this new and promising topic, which may lead to a new research field.

To make Bayesian inference on the considered model (1), we need to specify a prior distribution for unknown parameters or coefficients in approximating unknown nonparametric functions. A standard assumption for unknown parameters is some proper parametric distribution family such as normal distribution, gamma distribution, inverse gamma distribution, inverse Gaussian distribution, Wishart distribution, and Beta distribution in which their hyperparameters are user prespecified. Their misspecification or improper application may lead to unreasonable even misleading parameter estimation. Bayesian inference based on these assumptions did not utilize historical data and limits its popularity in that the usage of historical data may improve the efficiency of parameter estimation. To address this issue, some relaxed priors are considered, for example, see power prior, g-prior, normalized power prior [13], calibrated power prior, dynamic power prior, and the power-expected-posterior prior and the scale transformed power prior [14]. For a high-dimensional sparse parametric model, we can assume a spike and slab prior for the parameter, which can be hierarchically expressed as a mixture of a normal distribution and an exponential distribution.

Let ${\vartheta }_1$ be the set of unknown parameters associated with the distribution of ${\epsilon }_i$, and ${\vartheta }_2$ be the set of unknown parameters associated with the distributions of ${\epsilon }_i$ and ${\xi }_i$, ${\vartheta }_3$ be the set of unknown parameters associated with the distribution of ${\delta }_i$, and denote $\vartheta =\left\{{\vartheta }_1,{\vartheta }_2,{\vartheta }_3\right\}$. Let $\theta$ be a set of unknown parameters in $\left\{\beta ,\mathrm{\Lambda },\mathrm{\Pi },\mathrm{\Gamma },\vartheta \right\}$. Denote $Y=\left\{Y_i,i=1,\dots ,n\right\}$, $Y_{\mathrm{obs}}=\left\{Y_{i,\mathrm{obs}},i=1,\dots ,n\right\}$, $Y_{\mathrm{mis}}=\left\{Y_{i,\mathrm{mis}},i=1,\dots ,n\right\}$, $X=\left\{X_i,i=1,\dots ,n\right\}$, $F=\left\{{\omega }_i,i=1,\dots ,n\right\}$, where $Y_{i,\mathrm{mis}}$ and $Y_{i,\mathrm{obs}}$ are sub-vectors of $Y_i$ corresponding to missing and observed values, respectively. Let $D=\left\{Y_{\mathrm{obs}},X\right\}$ and $\delta =\{{\delta }_{\mathit{ij}},i=1,\dots ,n,j=1,\dots ,p$. If the marginal posterior distribution of $\theta$ given the dataset $D$ is

$$\begin{array}{l}\pi \left((|\theta |D\right)=\int \pi (\left.\theta ,F,Y_{\mathrm{mis}},\delta |D\right){\mathit{\text{dFdY}}}_{\mathrm{mis}}\mathit{d\delta }\\ =\frac{\int \pi \left(\theta \right)\pi \left(Y|\theta ,F,X\right)\pi \left(F|{\vartheta }_2\right)\pi \left(\delta |Y,{\vartheta }_3\right.{\mathit{\text{dFdY}}}_{\mathrm{mis}}\mathit{d\delta }}{\pi \left(Y_{\mathrm{obs}}|X\right)},\end{array}$$E2

its posterior mean (i.e., Bayesian estimate) can be evaluated by

$$\begin{array}{l}\widehat{\theta }=E\left.\left(\theta \right|D\right)=\int \mathit{\theta \pi }\left(\theta |D\right)\mathit{d\theta }=\frac{\int \mathit{\theta \pi }\left(\theta \right)\pi \left(Y|\theta ,F,X\right)\pi \left(F|{\vartheta }_2\right)\pi \left(\delta |Y,{\vartheta }_3\right){\mathit{\text{dFdY}}}_{\mathrm{mis}}\mathit{\text{dδdθ}}}{\pi \left(Y_{\mathrm{obs}}|X\right)},\end{array}$$E3

where $\pi \left(\theta \right)$ is the prior distribution of $\theta$, $\pi \left(Y|\theta ,F,X\right)$ is the probability density function of $Y$ given $\left(X,\beta ,\mathrm{\Lambda },{\vartheta }_1\right)$, i.e., the likelihood function of $\theta$ associated with the considered latent variable model, $\pi \left(F|{\vartheta }_2\right)={\mathrm{\Pi }}_{i=1}^n\pi \left({\eta }_i|{\xi }_i,\mathrm{\Pi },\mathrm{\Gamma },{\vartheta }_2\right)\pi \left({\xi }_i|{\vartheta }_2\right)$ is the probability density of $F$, $\pi \left(Y_{\mathrm{obs}}|X\right)=\int \pi \left(\theta \right)\pi \left(Y|\theta ,F,X\right)\pi \left(F|{\vartheta }_2\right)\pi \left(\delta |Y,{\vartheta }_3\right){\mathit{\text{dFdY}}}_{\mathrm{mis}}\mathit{\text{dδdθ}}$ is the marginal likelihood of $Y_{\mathrm{obs}}$ given $X$, and $\pi \left(\delta |Y,{\vartheta }_3\right)$ is the probability density of $\delta$ given $\left(Y,{\vartheta }_3\right)$.

From Eq. (3), it is easily seen that evaluating $\widehat{\theta }$ is almost impossible due to high-dimensional integral involved. To address this issue, the well-known Markov chain Monte Carlo (MCMC) algorithm is employed to approximate $\widehat{\theta }$ by sequentially drawing observations from the posterior distributions of components of $\theta$ and $\left(F,Y_{\mathrm{mis}},\delta \right)$ via the Gibbs sampler together with the Metropolis-Hastings algorithm. Denote $\pi \left(\theta |Y,F,X,\delta \right)$, $\pi \left(Y_{\mathrm{mis}}|Y_{\mathrm{obs}},F,X,\theta ,\delta \right)$, and $\pi \left(F|Y,\theta ,X\right)$ as the conditional distributions of $\theta$ given $(Y,F$, $X,\delta )$, $Y_{\mathrm{mis}}$ given $\left(Y_{\mathrm{obs}},F,X,\delta \right)$, and $F$ given $\left(Y,\theta ,X\right)$, respectively. The Gibbs sampler is implemented as follows. At the $t$th iteration of the Gibbs sampler with the current observations $\left\{{\theta }^{\left(t\right)},F^{\left(t\right)},Y_{\mathrm{mis}}^{\left(t\right)},{\delta }^{\left(t\right)}\right\}$ of $\left\{\theta ,F,Y_{\mathrm{mis}},\delta \right\}$, we sequentially draw i) ${\beta }^{\left(t+1\right)}$ from the conditional distribution $\pi \left(\beta |Y,F,X,\mathrm{\Lambda },{\vartheta }_1\right)$, ii) ${\mathrm{\Lambda }}^{\left(t+1\right)}$ from the conditional distribution $\pi \left(\mathrm{\Lambda }|Y,F,X,\beta ,{\vartheta }_1\right)$, iii) ${\mathrm{\Pi }}^{\left(t+1\right)}$ from the conditional distribution $\pi \left(\mathrm{\Pi }|F,\mathrm{\Gamma },{\vartheta }_2\right)$, iv) ${\mathrm{\Gamma }}^{\left(t+1\right)}$ from the conditional distribution $\pi \left(\mathrm{\Gamma }|F,\mathrm{\Pi },{\vartheta }_2\right)$, v) ${\vartheta }_1^{\left(t+1\right)}$ from the conditional distribution $\pi \left({\vartheta }_1|Y,X,\beta ,\mathrm{\Lambda },F\right)$, vi) ${\vartheta }_2^{\left(t+1\right)}$ from the conditional distribution $\pi \left({\vartheta }_2|F,\mathrm{\Pi },\mathrm{\Gamma }\right)$, vii) $Y_{\mathrm{mis}}^{\left(t+1\right)}$ from the conditional distribution $\pi \left(Y_{\mathrm{mis}}|Y_{\mathrm{obs}},\delta ,X,\beta ,\mathrm{\Lambda },F,{\vartheta }_1\right)$, and viii) $F^{\left(t+1\right)}$ from the conditional distribution $\pi \left(F|Y,X,\theta \right)$. The aforementioned conditional distributions may be some familiar distributions from which observations can be directly drawn. But in some cases, these conditional distributions may be some unfamiliar and rather complicated distributions from which observations can be indirectly drawn. In this case, some alternative approaches, for example, the Metropolis-Hastings algorithm, rejection sampling, acceptance-rejection sampling, importance sampling, hybrid-jump-based sampling, and reversible jump sampling, can be employed to sample observations from these complicated distributions. The convergence of the above introduced Gibbs algorithm can be monitored by the estimated potential scale reduction (EPSR) values associated with the parameters [15], which are evaluated continuously as the iterations run. The Gibbs sampler converges if the EPSR values of unknown parameters are less than 1.2. Also, we can assess the convergence of the Gibbs sampler by plotting several parallel sequences of observations drawn from different starting values of unknown parameters against iterations.

Let $\left\{{\theta }^{\left(m\right)},m=1,\dots ,M\right\}$, $\left\{F^{\left(m\right)},m=1,\dots ,M\right\}$ and $\left\{Y_{\mathrm{mis}}^{\left(m\right)},m=1,\dots ,M\right\}$ be $M$ observations sampled from their corresponding conditional distributions via the aforementioned Gibbs sampler after the Gibbs sampler algorithm converges, respectively. Bayesian estimates of $\theta$, $F$_, and $Y_{\mathrm{mis}}$ can be computed by

$$\widehat{\theta }=\frac{1}{M}\sum _{m=1}^M{\theta }^{\left(m\right)},\widehat{F}=\frac{1}{M}\sum _{m=1}^M{\theta }^{\left(m\right)},{\widehat{Y}}_{\mathrm{mis}}=\frac{1}{M}\sum _{m=1}^M{Y}_{\mathrm{mis}}^{\left(m\right)},$$E4

respectively. Their corresponding standard deviations can be computed with their corresponding sample covariance matrices of the observations. The details can refer to the literature [7]. The above argument on Bayesian inference is a classical method. However, for a high-dimensional parametric or nonparametric model, one needs some new approaches to solve the computing time and efficiency and stability of algorithm problem. In fact, when the dimension of covariate matrix is large and the sample size is relatively small, i.e., the well-known “large $p$ and small $n$” problem, the Gibbs sampler is computationally expensive and has poor stability.

To solve this issue for a high-dimensional regression model, there are some novel approaches developed for parameter/nonparametric function estimation in the Bayesian framework, for example, see Bayesian Lasso, Bayesian adaptive Lasso, Bayesian elastic net, and Bayesian $L_{1/2}$. These approaches can be utilized to estimate model parameters or nonparametric functions and are simultaneously used to select variables in a high-dimensional regression model, which have received considerable attention and extended to various models such as generalized linear models and linear mixed models. In particular, to reduce the computational cost, various variational Bayesian methods have been developed for various models in recent years. For example, see linear mixed models [16] and reference therein. However, there are a lot of unsolved problems. For example, for a complicated model, how to find the optimal variational densities for approximating complicated posterior distributions? How to extend/break the assumption of mean field that is a basic assumption in variational analysis? How to utilize other divergence criteria rather than Kullback–Leibler divergence to develop variational Bayesian theories?

3. Model comparison

Model comparison is widely used to select a plausible model to fit a given dataset among all the considered candidate models. Various methods have been developed to make model comparisons for many models such as linear/nonlinear regression models, structural equation models, multilevel models, machine learning models, and pattern recognition model in the Bayesian framework over the past years.

To select a better model among all the candidate models, we can adopt the well-known best subset selection methods such as Akaike information criterion (AIC), Bayesian information criterion (BIC), deviance information criterion (DIC), generalized information criterion (GIC), minimum description length (MDL), Hannan-Quinn information criterion (HIC), and log scoring criterion (also called the conditional predictive ordinate, i.e., CPO), which trade off a measure of model plausibility and a measure of model complexity. Also, the Bayes factor [17] has been developed to conduct Bayesian model comparison and is widely utilized to investigate the strength of the evidence in favor of one model among two candidate models. The Bayes factor for two competing models $H_0$ and $H_1$ is defined as follows:

$${\mathrm{B}}_{10}=\frac{\pi \left(Y_{\mathrm{obs}},\delta |H_1\right)}{\pi \left(Y_{\mathrm{obs}},\delta |H_0\right)},$$

where $\pi \left(Y_{\mathrm{obs}},\delta |H_k\right)=\int \pi \left(Y_{\mathrm{obs}},\delta |{\theta }_k\right)\pi \left({\theta }_k\right)d{\theta }_k$ is the marginal density of $H_k$ with parameter vectors ${\theta }_k$, and $\pi \left({\theta }_k\right)$ is the prior density of ${\theta }_k$ associated with model $H_k$ for $k=0,1$. In general, if the Bayes factor ${\mathrm{B}}_{10}>1$, the model $H_1$ is more plausible by the observed data than the model $H_0$, which leads to the following model comparison rule: ${\mathrm{B}}_{10}$‘s value lying in the intervals (3,10), (10,30), (30,100), and (100,$\infty$) yields moderate, strong, very strong, and extreme evidence in favor of model $H_1$, respectively. It is rather difficult to compute $\pi \left(Y_{\mathrm{obs}},\delta |{\theta }_k\right)$ due to the intractable high-dimensional integral involved, thus computing the Bayes factor ${\mathrm{B}}_{10}$ is challenging. Many methods have been proposed to compute marginal likelihoods $\pi \left(Y_{\mathrm{obs}},\delta |H_k\right)$ or Bayes factors [3]. For example, see importance sampling, path sampling, bridge sampling, Harmonic mean method, random weight importance sampling, sequential Monte Carlo method, and pareto-smoothed importance sampling leave-one-out cross-validation.

One serious defect of the Bayes factor for model comparison is that it is well defined for improper priors of ${\theta }_k$‘s and is sensitive to the selection of the hyperparameters in the priors. According to our experience, different priors together with different sampling methods lead to different values of the Bayes factor, i.e., different model comparison results. To this end, some modifications of the Bayes factor have been proposed, for instance, the partial Bayes factor, the intrinsic Bayes factor, and the fractional Bayes factor, which are subject to more or less arbitrary selection of training samples, weights for averaging training samples, and fractions, respectively. Also, some robust methods were developed to compute the sensitivity of the marginal likelihoods via the simulation-based methods, called the automated prior robustness method. Recently, some novel methods were proposed to deal with improper priors in computing the Bayes factor. For example, see machine learning method, i.e., first using a part of the dataset studied to train the Bayes factor/transform the improper prior into a proper prior and then utilizing the remainder of the dataset for model comparison, which provides a new idea for computing the Bayes factor. The robustness of model comparison is a challenging topic, which is worth further studying.