Sampling from Density power divergence-based Generalized posterior distribution
via Stochastic optimization

Let $\mathcal{X}\subset\mathbb{R}^{d}$ and $\Theta\subset\mathbb{R}^{p}$ denote a sample space and a parameter space, respectively. Suppose that $x_{1},\dots,x_{n}\in\mathcal{X}$ are independent and identically distributed from a distribution $F_{0}$ , and $\pi(\theta)$ is a prior probability density function. Let $f(\cdot\mid\theta)$ ($\theta\in\Theta$) denote the statistical model on $\mathcal{X}$. While Bayesian inference provides methods to quantify uncertainty in $\theta$, the true data-generating process $F_{0}$ is often unknown, and the assumed model may be misspecified.

To address this challenge, Bissiri et al., (2016) proposed a general framework for updating belief distributions through what is termed the generalized Bayesian inference. This framework enables quantifying the uncertainty of the general parameter $\theta$ which has the true value $\theta_{0}\coloneqq\mathop{\rm argmin}\limits_{\theta\in\Theta}\int\ell(\theta,x)dF_{0}(x)$, where $\ell(\theta,x)$ is a general loss function. The generalized posterior distribution takes the form:

where the loss function has an additive property: $\ell(\theta,x_{1:n})=\sum_{i=1}^{n}\ell(\theta,x_{i})$ (here and throughout, we use the notation $x_{1:n}$ to denote $(x_{1},\dots,x_{n})$), and $\omega>0$ is a loss scaling parameter that controls the balance between the prior belief and the information from the data. If we take $\omega=1$ and $\ell(\theta,x)=-\log f(x\mid\theta)$, then we get that $\pi(\theta\mid x)\propto\pi(\theta)f(x\mid\theta)$, which shows that the generalized posterior distribution includes the standard posterior distribution as a special case. The selection of the loss scaling parameter $\omega$ has been extensively studied, as it affects the posterior concentration and the validity of Bayesian inference (Grünwald and van Ommen,, 2017; Holmes and Walker,, 2017; Lyddon et al.,, 2019; Syring and Martin,, 2019; Wu and Martin,, 2023).

To achieve robust Bayesian inference against outliers, we can construct the loss function using robust divergence measures. For example, using DPD or $\gamma$-divergence (Fujisawa and Eguchi,, 2008) instead of the Kullback-Leibler divergence leads to inference procedures that are less sensitive to outliers while maintaining the theoretical guarantees of the generalized Bayesian inference (Ghosh and Basu,, 2016; Nakagawa and Hashimoto,, 2020). These robust divergences provide principled ways to deal with model misspecification and contaminated data while preserving the coherent updating property of Bayesian inference (Jewson et al.,, 2018).

The DPD offers a particularly attractive approach due to its robustness properties and theoretical guarantees. The DPD-based loss function (Basu et al.,, 1998) is defined for observations $x_{1:n}$ as:

where $\alpha>0$ is a tuning parameter and $F$ is the distribution corresponding to $f(\cdot\mid\theta)$. The DPD-based posterior (Ghosh and Basu,, 2016) provides a robust framework for updating prior distributions. The DPD-based posterior is defined as:

The DPD-based posterior inherits robust properties from the DPD loss function. Unlike traditional likelihood-based approaches, it is less sensitive to outliers and model misspecification. The robustness of the DPD-based posterior has been theoretically established through the analysis of influence functions (Ghosh and Basu,, 2016). The bounded influence function ensures protection against outliers and model misspecification, while maintaining good efficiency under the true model.

The tuning parameter $\alpha$ controls the trade-off between efficiency and robustness, where larger values of $\alpha$ provide stronger robustness at the cost of statistical efficiency under correctly specified models. The selection of this parameter in the Bayesian context has been extensively studied (Warwick and Jones,, 2005; Basak et al.,, 2021; Yonekura and Sugasawa,, 2023), providing both theoretical foundations and practical guidelines for choosing appropriate values based on the data characteristics and application context. These studies establish frameworks for parameter selection that maintain the balance between robust inference and computational tractability within the Bayesian paradigm.

The DPD-based posterior has robustness against outliers, which is a highly desirable property. However, the second term of the loss function (2) based on DPD contains an integral term (infinite sum) that cannot be written explicitly, except for specific distributions such as the normal distribution, making it difficult to calculate the posterior distribution using MCMC. While numerical integration can be used to approximate calculations in general probability models, numerical integration often fails for high-dimensional data (Niederreiter,, 1992). To address these computational challenges, we propose a new sampling method that avoids the need for explicit calculation or numerical approximation of the integral term. Before presenting our proposed approach, we first review in the next section the loss-likelihood bootstrap introduced by Lyddon et al., (2019), which will provide important theoretical foundations for our sampling algorithm. The loss-likelihood bootstrap offers a nonparametric way to approximate posterior distributions by resampling, but has limitations when applied directly to DPD-based inference. We will then build on these ideas to develop an efficient sampling scheme specifically designed for DPD posteriors with general parametric models.

The loss-likelihood bootstrap (LLB) (Lyddon et al.,, 2019) provides a computationally efficient method for generating approximate samples from generalized posterior distributions by leveraging Bayesian nonparametric principles. The key insight is to model uncertainty about the unknown true data-generating distribution through random weights on the empirical distribution. The algorithm proceeds as follows:

The asymptotic properties of LLB samples are characterized by Theorem 1 (Lyddon et al.,, 2019). This result shows that LLB samples exhibit proper uncertainty quantification asymptotically, with the sandwich covariance matrix $J(\theta_{0})^{-1}I(\theta_{0})J(\theta_{0})^{-1}$ providing robustness to model misspecification.

The LLB approach offers significant computational advantages compared to traditional MCMC methods. By generating independent samples and enabling trivial parallelization, it eliminates the need for convergence diagnostics while maintaining robustness to model misspecification. These properties make it particularly attractive for large-scale inference problems where computational efficiency is crucial.

Moreover, LLB plays a crucial role in calibrating the loss scale parameter for generalized posterior distributions. Under regularity conditions, both the generalized Bayesian posterior and LLB posterior distributions converge to normal distributions with identical means at the empirical risk minimizer $\hat{\theta}_{n}$ but different covariance structures: $w^{-1}J(\theta_{0})^{-1}$ for the generalized posterior and $J(\theta_{0})^{-1}I(\theta_{0})J(\theta_{0})^{-1}$ for the LLB posterior. By equating the Fisher information numbers, defined as $K(p)=\text{tr}(\Sigma^{-1})$ for a normal distribution with covariance $\Sigma$, they derive the calibration formula $w=\text{tr}\{J(\theta_{0})I(\theta_{0})^{-1}J(\theta_{0})^{\top}\}/\text{tr}\{J(\theta_{0})\}$. This calibration exhibits a remarkable theoretical property: when the true data-generating distribution belongs to the model family up to a scale factor (i.e., $f_{0}(x)=\exp\{-w_{0}\ell(\theta_{0},x)\}$ for some $\theta_{0}\in\Theta$ and $w_{0}>0$), it recovers the correct scale $w=w_{0}$. In practice, since $\theta_{0}$ is unknown, the matrices $I$ and $J$ are estimated empirically using the risk minimizer $\hat{\theta}_{n}$, resulting in a computationally tractable approach for calibrating the generalized posterior that ensures equivalent asymptotic information content with the LLB posterior. This theoretical result, established in (Lyddon et al.,, 2019), provides further justification for viewing LLB as an approximate sampling method from a well-calibrated generalized posterior distribution.

Recently, Okuno, (2024) made a significant breakthrough in minimizing the DPD for general parametric density models by introducing a stochastic optimization approach. While traditional DPD minimization required computing integral terms that could only be explicitly derived for specific distributions such as normal and exponential distributions, their stochastic optimization framework enabled DPD minimization for a much broader class of probability density models.

Inspired by this stochastic approach to minimizing DPD, we propose to reformulate the LLB algorithm (Lyddon et al.,, 2019) as a stochastic method. While LLB was traditionally formulated as an optimization problem based on weighted loss-likelihoods, introducing stochastic optimization techniques enables us to construct a more computationally efficient algorithm.

The objective function, $L_{w}(\theta)$, that we minimize in the LLB algorithm, is

We also calculate the gradient of the objective function:

where $u(\cdot\mid\theta)=\nabla\log f(\cdot\mid\theta)$. Then, we randomly generate independent samples $\xi^{(t)}=(z_{1}^{(t)},\dots,z_{m}^{(t)})$ from $F$, and define the stochastic gradient of step $t$ by

Taking the expectation of both sides of this equation, we have the following equation:

and confirm that (7) is an unbiased estimator. Let $v>0$ and assume the objective function $\mathcal{L}_{w}$ satisfies the following conditions: (i) $\mathcal{L}_{w}(\theta)$ is differentiable over $\Theta$ (ii) There exists $L>0$ such that for any $\theta_{1},\theta_{2}\in\Theta$, $|\nabla\mathcal{L}_{w}(\theta_{1})-\nabla\mathcal{L}_{w}(\theta_{2})|\leq L|\theta_{1}-\theta_{2}|$ (iii) $E_{\xi^{(t)}}[g(\theta^{(t)};\xi^{(t)})]=\nabla\mathcal{L}_{w}(\theta^{(t)})$ for all $t=1,\ldots,T$ (iv) $V_{\xi^{(t)}}[g(\theta^{(t)};\xi^{(t)})]\leq v$ for all $t=1,\ldots,T$. Under these conditions, as $T\to\infty$, the SGD algorithm converges to minimize the objective function $\mathcal{L}_{w}(\theta)$ in probability (Ghadimi and Lan,, 2013; Okuno,, 2024). Then, the parameter $\theta^{(t)}$ was iteratively updated by the SGD method:

where $\{\eta_{t}\}_{t=1}^{T}$ is the learning rate satisfying

Based on these theoretical foundations, we implement our approach through the SGD algorithm that offers two key advantages. First, it provides a unified framework that enables sampling from a broad class of parametric probability models, extending beyond the traditional limitations of explicit DPD computations. This generality makes the method applicable across diverse statistical scenarios without requiring model-specific derivations. Second, as demonstrated in subsequent simulation studies, the algorithm exhibits substantial computational efficiency even in high-dimensional settings, where traditional sampling methods often become intractable.

These dual characteristics - theoretical generality and computational tractability - establish our method as a significant contribution to modern statistical inference, particularly in settings requiring both model flexibility and efficient computation. We formalize this approach in Algorithm 2, which details the implementation of our stochastic gradient-based sampling procedure.

Building upon the stochastic gradient-based sampling framework introduced in the previous section, we now extend our methodology to encompass generalized linear models (GLMs). As in the previous section, our approach minimizes the objective function through SGD, and all theoretical properties, including the convergence conditions (i)-(iv), remain identical. This theoretical consistency allows us to leverage the same optimization framework while extending its application to the broader class of GLMs.

In GLMs, the random variables $y_{i}~{}(i=1,\dots,n)$ are independent and follow the general exponential family of distributions that have density:

where the canonical parameter $\theta_{i}$ is a measure of location depending on the fixed predictor $\mathbf{x}_{i}$ and $\phi$ is the scale parameter. Let $\mu_{i}=E[Y_{i}\mid\mathbf{x}_{i},\bm{\beta}]$, so that the GLM satisfies $h(\mu_{i})=\mathbf{x}_{i}^{T}\bm{\beta}$, where $h(\cdot)$ is a monotone and differentiable function.

Applying our stochastic gradient-based approach to GLMs, we consider Bayesian inference via DPD for the regression coefficients and the scale parameter. Let $\theta=(\beta,\phi)$ and let $\alpha>0$ be a hyperparameter, we can give the following loss function

Following the same principles as our general framework, we derive the stochastic gradient estimator for GLMs. The gradient of the objective function takes the form:

where $u(\cdot\mid\mathbf{x}_{i},\theta)=\nabla\log f(\cdot\mid\mathbf{x}_{i},\theta)$ and $F_{i}$ is the distribution corresponding to $f(y_{i}\mid\mathbf{x}_{i},\theta)$. We randomly generate independent samples $\xi_{i}^{(t)}=(z_{i1},\dots,z_{im})$ from the distribution $F_{i}$ and $\xi^{(t)}=(\xi_{1}^{(t)},\dots,\xi_{n}^{(t)})$. We can acquire the gradient of the objective function such that

Also, $g(\theta^{(t)}\mid\xi^{(t)})$ is an unbiased estimator from the following calculation.

This generalization advances DPD-based posterior sampling through stochastic gradient methods, particularly addressing models with intractable integral terms in the DPD. Unlike traditional approaches that rely on computationally expensive numerical integration, our method employs stochastic approximation, enabling efficient sampling for various GLMs including Poisson regression. Simulation studies demonstrate that this approach achieves superior computational efficiency and numerical stability while maintaining statistical accuracy. We formalize this method in Algorithm 3, which extends our framework to accommodate GLM-specific characteristics, particularly in gradient computation steps involving multiple samples per observation. The algorithm maintains the core structure of Algorithm 2 while incorporating these GLM-specific modifications.

We evaluated our proposed LLB algorithm, which replaces the objective function’s gradient with an unbiased stochastic gradient, through two comprehensive numerical experiments. For our evaluation dataset, we generated $1000$ samples from a normal distribution $\text{N}(0,1)$ and introduced outliers from $\text{N}(10,0.01)$ at a contamination ratio of $0.05$. We assumed a normal distribution with parameters $(\mu,\sigma)$ as our probability model. This setting is particularly advantageous as it allows us to compute the integral term of (2) explicitly as $(2\pi)^{-\frac{\beta}{2}}(1+\alpha)^{-\frac{3}{2}}\sigma^{-\alpha}$, enabling us to express the DPD-based posterior in closed form. This tractability facilitates the application of various sampling methods, including MCMC with the Metropolis-Hastings (MH) algorithm and the standard LLB algorithm, providing an ideal framework for comparative analysis.

Our first experiment validated that our stochastic gradient approach performs equivalently to the standard LLB algorithm. The results in Figure 1, based on 10000 posterior samples for each method, demonstrate that the performance of the LLB algorithm remains consistent regardless of the stochastic gradient implementation. In our proposed method, we obtained a posterior mean of $0.054$ with posterior variance $0.0013$ for $\mu$ and a posterior mean of $1.036$ with posterior variance $0.0007$ for $\sigma$. The standard LLB algorithm produced nearly identical results, with posterior mean $0.054$ and posterior variance $0.0013$ for $\mu$, and posterior mean $1.036$ and posterior variance $0.0007$ for $\sigma$.

Our second experiment confirmed the accuracy of our proposed method in sampling from the DPD-based posterior. Figure 2 presents overlaid histograms of samples ($10000$ samples for each method) obtained from the DPD-based posterior for $\mu$ and $\sigma$ using three methods: the MH algorithm, the Hamiltonian Monte Carlo (HMC) method implemented in Stan, and our proposed approach. Our method yielded posterior means of $-0.005$ and $0.988$ with posterior variances of $0.0012$ and $0.0007$ for $\mu$ and $\sigma$, respectively. These results closely aligned with those from the MH algorithm (posterior means of $-0.0057$ and $0.9893$, posterior variances of $0.0011$ and $0.0007$ for $\mu$ and $\sigma$) and the Stan implementation (posterior means of $-0.0054$ and $0.9899$, posterior variances of $0.0011$ and $0.0007$ for $\mu$ and $\sigma$). The close alignment of these distributions empirically validates that our method successfully samples from the target DPD-based posterior.

In high-dimensional settings, conventional sampling methods often face significant computational challenges. These challenges primarily manifest in increased computational time and reduced estimation accuracy. To validate the scalability and effectiveness of our proposed method in addressing these challenges, we conducted a comprehensive simulation study.

To evaluate the validity of our proposed method in high-dimensional settings, we designed our simulation framework using a contaminated multivariate normal model, where we generated $100$ samples from $N_{p}(\mathbf{0},\mathbb{I})$ with $5\%$ contamination from $N_{p}(10\mathbbm{1},0.01\mathbb{I})$. Here, $\mathbb{I}$ denotes the identity matrix and $\mathbbm{1}$ represents a vector of ones. For the numerical integration in the LLB algorithm, we performed the integration over a grid of $M^{p}$ points spanning the hypercube $[-D,D]^{p}$, where we set $M=10$ grid points per dimension and $D=2$ as the integration boundary. While this grid-based numerical integration approach provides high accuracy for 1-dimensional problems, the computational cost grows exponentially with dimension as it requires $M^{p}$ evaluations for $p$-dimensional integration. For the optimization in the LLB, we used standard gradient descent with a fixed learning rate $\eta=T^{-1}\sum_{t=1}^{T}\eta_{t}$. This setting allows us to systematically evaluate the effect of increasing dimensionality. We compared our proposed method against the LLB algorithm with numerical integration, both initialized with maximum likelihood estimates and configured to generate $1000$ posterior samples.

Table 1 compares the computational efficiency between our SGD approach and the LLB method using gradient descent with numerical integration (GD+NI) across varying dimensions. Through 1000 Monte Carlo simulations, we observed that both methods achieve similar mean squared error values for $p=2$ (SGD: 0.01281, GD+NI: 0.01463) and $p=3$ (SGD: 0.01341, GD+NI: 0.01476), indicating that the estimation accuracy is maintained regardless of the chosen method. The key difference lies in the computational time required to obtain 1000 samples. While SGD maintains consistently low computation times (approximately 3.6-4.0 seconds) regardless of dimension, GD+NI shows an exponential increase in computational cost as the dimension grows. The computational efficiency gain increases exponentially with dimension, reaching a 559-fold improvement at $p=7$ (3.63 seconds for SGD versus 2030.3 seconds for GD+NI). This substantial difference in computational time, while maintaining comparable estimation accuracy, clearly demonstrates the superior scalability of our stochastic approach in higher dimensions.

These empirical results not only support our theoretical findings but also demonstrate the practical viability of our method for high-dimensional robust Bayesian inference. The ability to maintain computational efficiency without sacrificing estimation accuracy makes our method particularly suitable for modern high-dimensional applications where traditional numerical integration approaches become computationally prohibitive.

To validate the practical utility of our method in scenarios where traditional approaches face computational challenges, we conducted extensive simulation studies using Poisson regression models. This choice of model is particularly important as it represents a common case where the integral term in (2) lacks a closed-form solution, thus requiring numerical approximation in conventional approaches.

For our simulation setup, we generated synthetic data comprising $n=300$ samples from a Poisson regression model, with the experiment repeated 100 times. The covariates $\mathbf{x}_{i}\in\mathbb{R}^{p}$ were drawn from a standard multivariate normal distribution $N_{p}(\mathbf{0},\mathbb{I})$, and the response variables were generated according to $y_{i}\sim\text{Po}(\exp(\beta_{0}+\mathbf{x}_{i}^{T}\bm{\beta}))$. The true parameter $(\beta_{0},\bm{\beta})$ was constructed by independently sampling each component from a uniform distribution over $[0,1/4]$, ensuring a realistic but challenging estimation scenario.

We compared our proposed stochastic optimization approach against the standard LLB algorithm that relies on finite sum approximation of the infinite sum. For the existing method, we employed the same gradient descent algorithm as described in the previous subsection. Following established practice in the literature (Kawashima and Fujisawa,, 2019), the LLB implementation used $M=n$ Monte Carlo samples for approximation. Both methods were initialized using maximum likelihood estimates and configured to generate $1000$ posterior samples, ensuring a fair comparison of their performance.

Table 2 presents a comprehensive comparison across varying dimensions ($p=2,5,10,\\ 20$) using three key performance metrics averaged over 100 simulation runs: mean squared error, coverage probability, and average length of 95% credible intervals. To evaluate performance, we first computed the posterior median $\hat{\beta}_{k}$ for each regression coefficient $\beta_{k}$ and assessed their accuracy using mean squared error, which is calculated as $p^{-1}\sum_{k=0}^{p}(\hat{\beta}_{k}-\beta_{k})^{2}$. Additionally, we constructed 95% credible intervals $CI_{k}$ for each $\beta_{k}$ and evaluated their effectiveness using two measures: the average length of the intervals, given by $p^{-1}\sum_{k=0}^{p}|CI_{k}|$, and the coverage probability, defined as $p^{-1}\sum_{k=0}^{p}I(\beta_{k}\in CI_{k})$. Our proposed method shows consistently superior performance across all dimensions, not only achieving lower MSE values but also maintaining coverage probabilities closer to the nominal 95% level. Furthermore, it provides more precise uncertainty quantification, as demonstrated by the shorter average interval lengths without compromising coverage.

These results convincingly demonstrate that our stochastic optimization approach delivers more accurate and reliable inference compared to traditional methods, particularly in settings where the model complexity precludes analytical evaluation of the objective function. The robust performance across increasing dimensions further highlights the scalability of our approach, making it particularly attractive for modern high-dimensional applications.