Adaptation of the tuning parameter in general Bayesian inference with robust divergence

We introduce a novel methodology for robust Bayesian estimation with robust divergence (e.g., density power divergence or γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}-divergence), indexed by tuning parameters. It is well known that the posterior density induced by robust divergence gives highly robust estimators against outliers if the tuning parameter is appropriately and carefully chosen. In a Bayesian framework, one way to find the optimal tuning parameter would be using evidence (marginal likelihood). However, we theoretically and numerically illustrate that evidence induced by the density power divergence does not work to select the optimal tuning parameter since robust divergence is not regarded as a statistical model. To overcome the problems, we treat the exponential of robust divergence as an unnormalisable statistical model, and we estimate the tuning parameter by minimising the Hyvarinen score. We also provide adaptive computational methods based on sequential Monte Carlo samplers, enabling us to obtain the optimal tuning parameter and samples from posterior distributions simultaneously. The empirical performance of the proposed method through simulations and an application to real data are also provided.


Introduction
One well-known way to deal with outliers and model misspecification when conducting inference is to use robust divergences.Since the pioneering work of Basu et al. (1998) that proposed density power divergence as an extension of the standard likelihood, some variants of the divergence (e.g.Fujisawa and Eguchi, 2008;Cichocki et al., 2011;Ghosh et al., 2017) and various statistical methods using robust divergences have been developed.Many robust divergences are indexed by a single tuning parameter that controls the robustness against outliers.If the tuning parameter is set to a smaller value than necessary, the resulting estimator may still be affected by outliers.On the other hand, using an unnecessarily large value for the tuning parameter leads to a loss of statistical efficiency (Basu et al., 1998).
Despite the success of the theoretical analysis of properties of statistical methods based on robust divergences, how to adaptively estimate the tuning parameter from the data has often been ignored, with a few exceptions (Warwick and Jones, 2005;Basak et al., 2021) that propose frequentist methods to select the optimal value via asymptotic mean squared errors.
There is a growing body of literature on Bayesian approaches using robust divergence.For example, general theory has been considered by (Ghosh and Basu, 2016;Jewson et al., 2018;Nakagawa and Hashimoto, 2020) and some specific applications to statistical models such as linear regression (Hashimoto and Sugasawa, 2020), change point detection (Knoblauch et al., 2018) and Bayesian filtering (Boustati et al., 2020).Nevertheless, a reasonable estimation strategy for the tuning parameter has not been carefully discussed.A natural consideration to find the best tuning parameter in the context of Bayesian statistics will be the use of model evidence or marginal likelihood.However, as we shall illustrate later, evidence is not useful for choosing the tuning parameter since the exponential of robust divergence cannot be directly interpreted as a normalised statistical model.
In this paper, this issue is addressed by taking advantage of ideas from statistical theories for unnormalised statistical models (Hyvärinen, 2005) introducing the Hyvarinen score (H score), which is a finite version of Fisher divergence.Based upon the idea of Hyvärinen (2005), Dawid et al. (2015); Shao et al. (2019) consider unnormalisable marginal likelihoods, with particular attention to model selection, where such unnormalisability is driven by the improper prior.Our main idea is to regard the exponential of robust divergence as unnormalisable models and employ a posterior distribution driven by the H-score, inspired by Dawid et al. (2015); Shao et al. (2019), as an objective function of γ.Since our objective function cannot be computed analytically in general, we then take advantage of sequential Monte Carlo (SMC) samplers (Del Moral et al., 2006) within a Robbins-Monro stochastic gradient framework to approximate the objective function, to estimate the optimal tuning parameter and to obtain samples from posterior distributions.Therefore, our work can be understood as an attempt to fill the current gap between the existing theory of such robust estimation and their practical implementation.Our proposed method has the following advantages over existing studies.i) Unlike existing methods (Warwick and Jones, 2005;Basak et al., 2021), our proposed method does not require a pilot plot.To optimise a tuning parameter, it is necessary to determine a certain value as a pilot estimate.Thus, the estimates may be strongly dependent on the pilot estimate.In addition, such methods often estimate excessively large values of the tuning parameters, given the proportion of outliers in the data.In contrast, our algorithm is stable and statistically efficient since that does not require a pilot estimate.
ii) Proposed methods in (Warwick and Jones, 2005;Basak et al., 2021) require an analytical representation of the asymptotic variance that cannot be obtained in general.
Compared to such methods, our proposed method does not require such expression, and therefore our method can be applied to rather complex statistical models, which seem difficult to be handled by the previous methods.
iii) We take advantage of SMC samplers (Del Moral et al., 2006) in the generalised Bayesian framework (Bissiri et al., 2016) with a gradient-ascent approach to perform parameter inference for the tuning parameter and posterior sampling simultaneously.This is a unique favourable algorithmic characteristic compared to methods that estimate parameters by fixing tuning parameters or by estimating tuning parameters once and then estimating other objects of interest.
Recently, Jewson and Rossell (2021) 2021) is that we use a natural form of general posterior based on robust divergence, which is widely adopted in the literature (e.g.Ghosh and Basu, 2016;Jewson et al., 2018), while the form of the posterior distribution in Jewson and Rossell (2021) is different from ours.As we mentioned, the main contribution of their work is the construction of model selection criteria of the BIC type through the Laplace approximation and the proof of their consistency.On the other hand, our research is about the estimation of tuning parameters and the inference of the posterior distribution, and the main objective is to propose an objective function and a computational method considered suitable for this purpose.In the framework of Generalised Bayesian, some methods that use variational Bayesian inference (Knoblauch et al., 2019;Frazier et al., 2021) have also been proposed.However, instead of computational speed, the approximate distribution obtained does not match the target distribution in the limit, and the method of estimating the tuning parameters is unclear.There is also a need to make some natural but somewhat stronger assumptions than our proposed method, such as that the target distribution is an exponential family.
The rest of the paper is organised as follows.In Section 2, we first set up the framework and then show theoretically and numerically that evidence induced by density power divergence to select the tuning parameter.Instead, we propose to estimate it based on the H-score (Hyvärinen, 2005;Dawid et al., 2015) and characterise its asymptotic behaviour.
As mentioned earlier, our method involves functions for which it is difficult to obtain an analytic representation.Therefore, we develop an adaptive and efficient Markov chain Monte Carlo (MCMC) algorithm based on SMC samplers (Del Moral et al., 2006) in Section 3. Numerical applications of the developed methodology are provided in Section 4, then conclusions and directions for future research are provided in Section 5.
2 Bayesian Inference with Robust Divergence ∼ G where G denotes the true distribution or the data-generating process.Also, assume that we have a statistical model {f θ : θ ∈ Θ} where Θ ⊆ R d for some d ≥ 1.We then write y 1:n to denote (y 1 , . . ., y n ) and let g denote the density of G with respect to dy.To make robust Bayesian inferences for θ, we use a potential function based on robust divergence instead of the standard likelihood function.
Here, we simply consider the density power divergence (Basu et al., 1998) but other ones with tuning parameters, such as γ-divergence (Fujisawa and Eguchi, 2008;Nakagawa and Hashimoto, 2020), α-divergence (Cichocki and Amari, 2010), Hölder divergence (Nielsen et al., 2017), also can be used within our framework.Given a prior density π(θ) with respect to dθ, we can define the corresponding posterior density where γ ∈ (0, ∞], p γ (y 1:n ) := Θ L γ (y 1:n ; θ)π(θ)dθ and log L γ (y 1:n ; θ) := (2.2) Note that p γ (y 1:n ) is generally referred to as evidence.In many scenarios, robustified posterior densities such as (2.1) give much more accurate and stable inference against outliers and theoretical properties of the posterior have been investigated in (Ghosh and Basu, 2016;Nakagawa and Hashimoto, 2020).However, its performance depends critically on the choice of the tuning parameter γ in (2.1) (e.g.Ghosh and Basu, 2016), which motivate us to find "best" γ to make inference successful.Notice that (2.1) can be seen as a special case of general Bayesian updating (Bissiri et al., 2016) with weight setting 1.As noted in Jewson et al. (2018), the density power divergence does not have any arbitrariness in the scale as a loss function, and one can set ω = 1.Under the general framework of Bayesian updating, Corollary 1 of Fong and Holmes (2020) implies that evidence is still the unique coherent marginal score for Bayesian inference.Thus, from the viewpoint of Bayesian statistics, it appears to be natural to find the best γ based on evidence, but its property of it is unclear since L γ (y 1:n ; θ) is not a probability density of y 1:n .Furthermore, the tuning parameter γ cannot be interpreted as "model parameter" in this case.The following example highlights the problem of using p γ (y 1:n ) to find the best γ.

Failure of model evidence: motivating example
To see why evidence is not useful for estimating γ, we start with the following proposition for a rescaled log L γ (y i ; θ).
Proof.See Appendix A.
Since the term − 1 γ + 1 is eliminated when considering the posterior distribution (2.1), this rescaling is a non-essential modification in the method we shall propose later.The meaning of the rescaling is to ensure that log L R γ (y i ; θ) converges to the log-likelihood as γ → 0, so that log L R γ (y i ; θ) can be regarded as a natural extension of the log-likelihood.
The important point here is that Proposition 1 implies that there are theoretically at least some situations where evidence is increasing monotonically for γ.Indeed, the following numerical example vividly illustrates such a situation.This implies that one cannot estimate γ using p γ (y 1:n ).A similar phenomenon is also discussed in Jewson and Rossell (2021).

Estimation using H-score
To overcome the illustrated problem, we first treat L γ (y 1:n ; θ) as an unnormaliseable statistical model motivated by Hyvärinen (2005).Note that even with an unnormalisable model, the update in (2.1) can be considered as a valid belief update according to Bissiri et al. (2016).It should be noted that when log L γ (y 1:n ; θ) is the density power divergence of the form (2.2), the normalising constant may not exist.For example, when f θ is a normal distribution, log L γ (y i ; θ) converges to a constant value under |y i | → ∞, so the integral of L γ (y 1:n ; θ) with respect to y 1:n diverges.Recently, Jewson and Rossell (2021) has pointed out that the role of such unnormalisable models can be recognised in terms of relative probability.
For d y dimensional observations y and twice differentiable density p(•), Hyvärinen (2005) defines the H-score as We then select the optimal γ with the smallest leave-one-out H-score, defined as Shao et al. (2019) adoptes the H-score to define prequential score for state space models, and the criteria (2.3) can be seen as prequential score under i.i.d.settings.As shown in Appendix B, under the assumptions stated in Shao et al. (2019), the leave-one-out H-score (2.3) can be rewritten as where the expectation is with respect to the robustified posterior distribution (2.1).Then, we can estimate γ as follows γ := arg min γ H n (γ). (2.5) As we shall discuss later, it can be shown that, under some conditions, fore, H n (γ) can be considered as an empirical approximation of the Fisher divergence J (γ) for the marginal distribution based on unnormaliseable models defined by robust diver- gence.An important point here is that the estimation by the H-score is independent of the normalisation constant.The following proposition is theoretical justification of selecting γ via (2.5).
Proof.See Appendix C.
Remark 1.As we mentioned, the prequential version of H n (γ) is also called the H-score in the context of Bayesian model selection (Shao et al., 2019;Dawid et al., 2015).The main advantage of using the H-score in this context is that it will provide a consistent and coherent model selection criterion.Jewson and Rossell (2021) proposes a consistent model selection criterion that is similarly based on H-scores but with batch estimation.
Although the prequential method is coherent, this comes with a very high computational cost, for every model, one must do posterior inference on all permutations of the data and increasing sample sizes.Here, we use a batch estimation approach to estimate γ, which avoids high computational costs.We also want to emphasise that, as we shall study later, such a batch approach will give rise to natural and efficient algorithms to estimate γ and posterior sampling.
Under the density power divergence (2.2), the first and second order derivatives of log L γ (y i ; θ) are given by These expressions do not include the integral term f θ (x) 1+γ dx, which makes the calculation of H n (γ) much more straightforward in practice since the integral term often is a form of a complicated expression.

Numerical illustration of the H-score under normal distribution
We consider the same problem in the example in Subsection 2.2.For a normal distribution N (µ, σ 2 ), the derivatives are as follows where φ(•; µ, σ 2 ) is the density function of N (µ, σ 2 ).We calculated H n (γ) in (2.4) for each γ, where posterior expectations were approximated by 2000 posterior samples of µ.
The data were simulated in the same way as in Subsection 2.2.The results are shown in Figure 2.2 when τ = 10 (blue lines) and 30 (red lines).Our experiment shows numerically that H n (γ) has a local minimum.Furthermore, it can be seen that in regions where γ is small, the posterior mean is relatively heavily influenced by outliers.In contrast, the posterior mean settles to a constant value in regions where γ is greater than the value that minimises H n (γ).Uncertainty in the sense of CI becomes greater.This result would suggest that statistical inefficiencies occur in regions where γ is larger than is necessary.

Sequential Monte Carlo Samplers
A natural way to obtain γ in (2.5) will be to use Robbins-Monro-type recursion.
where t κ t = ∞, t κ 2 t < ∞ but, in general, posterior sampling based on the Monte Carlo approximation will be required to evaluate ∇ γ H t (γ t ).That is, we need to construct an estimator of ∇ γ H t (γ t ).To do so, we first treat γ t in (3.1) as the positive sequence such gives rise to the following tempering-like distributions on a common measurable space, say (Θ, B(Θ)) The meaning of tempering-like here is not the usual tempering.It means that a family of distributions is constructed in the same measurable space by a sequence of updated γ t , gradually approaching the target distribution in the sense of an optimised γ, say γ .Using this, we can define a sequence of distributions defined on product spaces (Θ t , B(Θ t ) := where L k is a transition kernel from Θ k+1 to Θ k .Notice that Πt (θ 0:t | y 1:n ) admits marginally Then it is given by by construction, as the Radon-Nikodym derivative between them, one can derive unnormalised incremental weights as follows log w j t := log Given the time steps T > 0 and initial γ 0 > 0, we update γ t via (3.5) and iterate the SMC sampler T times.Our method can be algorithmically summarised as follows.
iii) Calculate w j t in (3.4) and set Remark 2. Since {w j t } are independent of {θ i t } but dependent of {θ i t−1 }, the particles {θ Theoretical guarantees of convergence of the Robbins-Monro algorithm usually require et al. (2016) shows that ∇ γ Ĥt (γ t ) is still a (weakly) consistent estimator, but not an unbiased estimator.Such unbiased estimation may be possible by using the recently developed MCMC with couplings in Algorithm 1, see Middleton et al. (2019); Jacob et al. (2020) for details.Instead of discussing convergence through such unbiased estimation, we shall discuss convergence through numerical experiments in the following sections.
We end this section by noting several advantages of the proposed method.First, Algorithm 1 enables us to estimate the tuning parameter and obtain posterior sampling simultaneously.This is a notable difference from existing methods, such as running MCMC or the EM algorithm with a fixed tuning parameter, for example, Fujisawa and Eguchi (2008); Ghosh and Basu (2016).We believe that it may be emphasised that by setting up a well-defined target function and using the stochastic gradient framework-based SMC samplers, it is possible to avoid the two-stage estimation that many previous studies have done in this context.We also emphasise that our proposed method has two notable advantages over existing methods: it does not require pilot plots, and it does not require an expression of the asymptotic variance of the model.Next, recall that as γ ↓ 0, L γ (y 1:n ; θ) converges to Kullback-Leibler divergence.Let γ be the value of converged γ in Algorithm 1. Then Algorithm 1 may be producing an approximated bridge between (multiplied by a prior distribution) Kullback-Leibler divergence and the target distribution L γ (y 1:n ; θ).
Therefore, sampling from such tempering-like distributions induced by the density power divergence (3.2) could provide a beneficial tempering effect and a potential reduction in computational complexity, particularly when d is large (Neal, 2001).Finally, Algorithm 1 will give rise to a natural way to construct adaptive MCMC kernels.Suppose that we use Metropolis-Hastings kernels based on a Gaussian random walk.Notice that before MCMC step in Algorithm 1, we have available.Also Ghosh and Basu (2016, Theorem 1) shows that Π t−1 (θ | y 1:n ) can be approximated by the Gaussian distribution.These will lead us to set M t (θ ∼ N (0, 2.38d −1/2 Σt−1 ), for instance.We note that this proposal is also can be considered as a consequence of the results from the optimal scaling analysis for random walk
ADAM is nowadays a standard and very effective addition to the type of recursive inference algorithms we are considering here, even more so as for increasing dimension of unknown parameters.See the above references for more motivation and details.

Simulation Studies
We here demonstrate the numerical performance of the proposed method.Throughout this study, we consider Gaussian models {f θ : θ ∈ Θ} = N (µ, σ 2 ) where both µ and σ 2 are unknown parameters.We first simulated data {y i } 100 i=1 i.i.d.
∼ N (1, 1) and then randomly replaced τ % of {y i } 100 i=1 by y i + 5.This setting is commonly referred to as the M-open world (Bernardo and Smith, 2009) in the sense that there is no θ such that g = f θ .

Experiment 1: Convergence property
We investigate convergence behaviour of Algorithm 1.In this study, we set (N, T, γ 0 ) = (2, 000, 500, 0.1) with 50 MCMC steps in Algorithm 1 to estimate γ.The results are shown in Figure 4.1.As can be seen in the figure, our proposed method converges stably to the true value after about 100 iterations.Here, the true value was obtained by first approximating the H-score with MCMC in the same way as before and then using a grid search to find the γ that minimises it.Experiment 2: Comparison with methods using fixed values of γ We next compare the performance of our proposed method with that of a non-adaptive method using a fixed value of γ.We set τ ∈ {0, 10, 20, 30}, and for each case we computed the posterior distribution of µ using Algorithm 1 and the vanilla version of MCMC.
For Algorithm 1, the tuning parameters were set to (N, T, γ 0 ) = (2, 000, 300, 0.1) with 50 MCMC steps, and for vanilla MCMC, γ was set to γ ∈ {0.1, 0.3, 0.5, 0.7, 0.9} in advance of the estimation with 100,000 MCMC steps.We used Metropolis-Hastings kernels based on a Gaussian random walk with N (0, 0.4) for the two cases under the uniform prior.Using the posterior samples obtained from Algorithm 1 and non-adaptive methods, we computed the posterior mean and 95% credible interval of µ.We ran 100 Monte Carlo experiments to calculate their (empirical) mean square error (MSE) and average 95% credible interval (ACI).The MSE is computed against the target value of 1, and the value is multiplied by 100.The results are given in Table 1.Although it is a simple example, the results summarised in the table clearly show that the accuracy of the inference is improved by estimating γ from the data rather than simply fixing it in terms of MSE.In fact, the best γ among the five choices depends on the underlying contamination ratio that we do not know in practice.Hence, it is difficult to determine a suitable value of γ simply by looking at the data, while our method can automatically tune the value of γ from the data.It should also be noted that the importance of adaptive tuning of γ is reflected in the results of not only MSE but also ACI; that is, the interval length obtained from Algorithm 1 is narrow compared with the non-adaptive methods in all the four scenarios.We next compare the proposed method with the H-posterior proposed by Jewson and Rossell (2021) (denoted by JR hereafter), where the posterior of the model parameters (µ, σ 2 ) as well as γ can be obtained.To apply the JR method, we generated 1000 posterior samples after discarding the first 500 samples.We evaluate the performance of the inference of µ and σ by MSE (multiplied by 100), coverage probability (CP) and average length (AL) of 95% credible intervals.The results are shown in Table Table 2, where the average estimates of γ are also shown.Although both methods provide similar estimates of γ, the accuracy of point estimation of JR is slightly better than that of Algorithm 1.However, it is observed that JR tends to produce a short coverage length so that the CP of the JR method is much smaller than the nominal level 95%.Accordingly, the average length is much smaller than those by Algorithm 1.This means that a direct application of the Hposterior by Jewson and Rossell (2021)

Experiment 4: Comparison with tempering
Following Nakagawa and Hashimoto (2020), we compare the robustness to outliers for the two generalised posterior distributions.The first distribution is constructed in the same way as before, while the other is constructed using tempering.We specify a tempered posterior the sequence, we divided the interval [0, 1] into 500 equal parts.We applied Algorithm 1 and the SMC sampler with the tempered posterior to test the robustness of the proposed method to data sets containing outliers.We used N = 2, 000 particles, 50 MCMC steps for both methods, and set (T, γ 0 ) = (500, 0.1) for Algorithm 1.The prior and MCMC kernels were set as the previous experiment, and the density estimation results obtained from the estimation results are summarised in Figure 4.2.The red line represents the posterior density estimate for µ when the data do not contain any outliers, and the blue line represents it when the data contain outliers.It is clear from the estimation results that our proposed method is robust even when the dataset contains outliers, while the SMC sampler with the tempered posterior is greatly affected by outliers, and the estimated posterior distributions are completely separated as a result.

Newcomb data
We apply our methodology to Simon Newcomb's measurements of the speed of light data, motivated by applications in Stigler (1977); Basu et al. (1998); Basak et al. (2021).The data can be obtained from Andrew Gelman's webpage: http://www.stat.columbia.edu/~gelman/book/data/light.asc.The sample size of the data set is 66 and contains two outliers, -44 and -2, illustrated in Figure 4.3.We fitted a Gaussian distribution model {f θ : θ ∈ Θ} = N (µ, σ 2 ) to the data and used Algorithm 1 to obtain the posterior distribution of the parameters (µ, σ).The tuning parameters (N, T, γ 0 ) in Algorithm 1 were set to (2, 000, 300, 0.1) with 50 MCMC iterations.The MCMC kernel was constructed as in the previous examples, and results are given in Figure 4.4.The existing study (Basak et al., 2021) reported γ = 0.23 for the same data set, which is very high compared to our estimate result of γ = 0.0855.Since the method proposed in Basak et al. (2021) requires a pilot plot and the estimation results depend significantly on it, we believe our estimation results are more reasonable.In fact, it is unlikely that we will have to use a value of γ = 0.23 for a data set that contains only two outliers.As shown in Basu et al. (1998), the parameter estimates are almost the same when γ = 0.0855 and when γ = 0.23.However, from the point of view of statistical efficiency, it would be preferable to adopt the lower value of γ = 0.0855 if the estimates were the same.To confirm this, 100 bootstrap resamplings were performed on the data, and the posterior bootstrap mean of each parameter and the variance was calculated, reported in Table 3.For each re-sampled data, we compared the results when the posterior distribution was calculated while estimating γ with our method and when the posterior distribution was calculated using MCMC after estimating and fixing it with the method proposed in Basak et al. (2021).The numerical experiments show that although the means estimated parameters agree between the two methods, the variances are much smaller for our method, suggesting that overestimation of γ leads to statistical inefficiency.3: Mean and variance of the posterior means of the parameters from 100 bootstrap re-samplings.The first row shows the result when using the proposed method, and the second one shows when using the method studied in Basak et al. (2021)

Hertzsprung-Russell Star Cluster Data
We use our methodology to perform linear regression models with normal errors, that is Basak et al. (2021), we fitted the regression model to the Hertzsprung-Russell star cluster data (Rousseeuw and Leroy, 2005), without constants.
The data set contains 47 observations on the logarithm of the effective temperature at the surface of the CYG OB1 star cluster (Te, covariates {x i }) and the logarithm of its light intensity (L/L0, explained variables {y i }).The data can be obtained from https: //rdrr.io/cran/robustbase/man/starsCYG.html, and shown at Figure 4.5.The tuning parameters (N, T, γ 0 ) in Algorithm 1 were set to (2, 000, 300, 0.1) with 50 MCMC iterations, and we used the uniform prior for (β, σ).The MCMC kernel was constructed as in the previous examples, and results are given in Figure 4.6.The corresponding OLS estimates were ( β, σ) = (1.1559,0.7219).Whilst we obtained γ = 0.1165, Basak et al. (2021) reported γ = 0.76 for the same data set.Our numerical experiments and previous studies (Ghosh and Basu, 2016;Nakagawa and Hashimoto, 2020) will suggest that as the proportion of outliers in the data increases, the value of γ also tends to increase.Thus, such a large value of γ is not reasonable considering the proportion of outliers in the data (only four samples in the lower right part in Figure 4.5), suggesting the superiority of our proposed method.Indeed, to confirm the suggested statistical inefficiency, the same experiments as in the previous section were carried out, and the results are summarised in Table 4.Although the results are not as striking as in the previous Newcomb data example, it would be possible to confirm that statistical inefficiencies occur in the estimation by the proposed method in   4: Mean and variance of the posterior means of the parameters from 100 bootstrap re-samplings.The first row shows the result when using the proposed method, and the second one shows when using the method studied in Basak et al. (2021) 5 Concluding remarks Our proposed method performs reasonably well and provides one of the few options, as far as we know, for routine robust Bayesian inference.To the best of our knowledge, this is the first attempt to propose both a theory and a computational algorithm to estimate the tuning parameters from data and to allow robust Bayesian estimation.We have shown numerically that a more efficient Bayesian estimation can be achieved by estimating the tuning parameter γ from the data.We have proposed an efficient sampling method using SMC samplers considering the sequence of γ as the temperature.Compared to existing studies (Warwick and Jones, 2005;Basak et al., 2021), our method has specificity and usefulness in that we can estimate the tuning parameters and sample from the posterior distribution simultaneously, and pilot plots and the asymptotic variance formula are not necessary.In this paper, we have focused in particular on the case of the density power divergence, but we want to stress that our method is general enough in the sense that it can be applied to the Bayesian estimation of other robust divergence-induced models.Furthermore, our framework opens up a number of routes for future research and insight, including those described below.
i) As we have noted, the integral term f θ (x) 1+γ dx is eliminated in the H-score, while the computation of the posterior distribution is computationally expensive, so it may be better to consider the H-posterior studied in Jewson and Rossell (2021) in this respect.However, the robustness of the H-posterior has not yet been studied, and it would therefore be interesting in the future to investigate this point in more detail using the influence function.
ii) Another direction of investigation involves the construction of an efficient MCMC kernel for posterior distributions derived from robust divergence such as (2.1).To make good inferences from data containing outliers, the posterior distribution induced by robust divergence is a model with a more or less heavy-tailed.As studied in Kamatani (2018), many standard MCMC algorithms are known to perform poorly in such cases, especially in higher dimensions.Therefore, studying MCMC algorithms within Algorithm 1 tailored to the posterior distribution induced by robust divergence would allow for more efficient Bayesian robust estimation.
iii) In this study, we have not focused on time series data, but, as Shao et al. (2019) shows, the H-score can also be defined for models that are not independent, for example, state-space models.In fact, Boustati et al. (2020) proposes a method for Bayesian filtering of state-space models using robust divergence, but the tuning parameters need to be estimated before filtering, and in this sense, it is not online filtering.Algorithm 1 does batch estimation, but we believe that extending it to online estimation would allow robust filtering of the state-space model while estimating the tuning parameters online from the data.iv) This study has focused on estimation and computational methods proposed in the generalised Bayesian framework, particularly using robust divergence.On the other hand, methods using the Maximum Mean Discrepancy (Chérief-Abdellatif and Alquier, 2020) and Kernel Stein Discrepancy (Matsubara et al., 2022) have also been proposed in recent years in the same generalised Bayesian framework, although not with the motivation of dealing with outliers.Both require adjustment of the hyperparameters of the kernel used, and it may be possible to estimate them using our proposed method and compare their performance.We have avoided comparing these potential alternative approaches because we believe this would obscure the main messages we have tried to convey within the numerical results section.Such a detailed numerical study can be the subject of future work.Shao et al. (2019), the following identify always holds where the expectation is taken with respect to the posterior distribution of θ given y.Using the above identity with p(y) = p(y i |y −i ), p(y|θ) = p(y i |θ, y −i ) and Π(θ) = Π(θ|y −i ), where Notice that, since we have assumed i.i.d.observations, we have p(y i | y i−1 , θ) = p(y i | θ).

C Proof of Proposition 2
The proof here is essentially the same as , see Dawid et al. (2015); Shao et al. (2019).Then under the assumptions stated in supplement of Shao et al. (2019), the variance term will converge at 0 w.p.1.Let (B, • ) be the space of continuous real functions on the compact set of γ equipped with the sup-norm.
Then, under the same assumptions, ] may take values in this space.Then strong law of large numbers on a separable Banach space (Azlarov and Volodin, 1982;Beskos et al., 2009), may be applied to Combining this result with (s10) in supplement of Shao et al. (2019) and integration by parts (Hyvärinen, 2005;Dawid et al., 2015) would yield lim n sup γ 1 n H n (γ) = J (γ) w.p.1.The result follows under the assumption for identification such that J (γ) is only maximised at γ , see A1. in Jewson and Rossell (2021) for instance.

D Derivatives of the H-score
In the following argument, we assume the exchangeability of integral and derivative without any remarks.For simplicity, we consider a univariate case, namely d y = 1.We define D γ (y 1:n ; θ) := log L γ (y 1:n ; θ) and D γ (y i ; θ) := log L γ (y i ; θ).The derivative of the H-score with respect to γ is expressed as

D.2 Normal distribution case
When y i ∼ N (µ, σ 2 ), the corresponding density power divergence is The detailed expressions of quantities appeared in the derivative of the H-score in (D.4) are obtained as follows: where w i = φ(y i ; µ, σ 2 ) γ .When y i ∼ N (x i β, σ 2 ), the derivative of the H score in the model is obtained by replacing µ with x i β.

Figure 2
Figure 2.1: Estimated p γ (y 1:n ).The Y-axis represents the value of p γ (y 1:n ), and the X-axis represents the value of γ.

Figure 2
Figure 2.2: The top left-hand plots the values of H n (γ) on the Y-axis, and values of γ on the X-axis when τ = 10.The minimum value of H n (γ) was obtained when γ = 0.1874.The top right-hand figure plots the sample mean values of the posterior mean of µ on the Y-axis and the corresponding values of γ on the X-axis under the same setting.The vertical line represents the value of γ that minimises H n (γ), and the thin ribbon line represents the 95% credible interval.The bottom left-hand and right-hand figures represent the same figures when τ = 30 respectively.In this case, he minimum value of H n (γ) was obtained when γ = 0.3311.
after resampling in Algorithm 1.In addition, Algorithm 1 uses a simple multinomial resampling applied at each step.The variability of the Monte Carlo estimates can be further reduced by incorporating dynamic resampling via the use of effective sample size.See Del Moral et al. (2006); Dai et al. (2020) for details.

Figure 4 . 1 :
Figure 4.1: Gaussian models experiment: Trajectories from execution of Algorithm 1.We used N = 2, 000 particles with 50 MCMC iterations and initial value γ 0 = 0.1.The left panel shows results when τ = 5 and the right one shows when τ = 10.The horizontal dashed lines in the plots show the true parameter γ = 0.2339 for the left and γ = 0.3638 for the right.The blue lines show the trajectory of γ estimated by Algorithm 1.

Figure 4 . 2 :
Figure 4.2: The density estimation of µ estimated by Algorithm 1 (left) and SMC sampler with the tempered posterior (right).The blue line shows when τ = 20 (contains 20% outliers) and the red one shows when τ = 0 (contains 0% outliers) in both panels.

Figure 4
Figure 4.3: The histogram of Simon Newcomb's measurements of the speed of light data.

Table 1 :
Jewson and Rossell (2021)rors (MSE) and average 95% credible intervals (ACI) of Algorithm 1 and the non-adaptive method (the vanilla version of MCMC with fixed γ), based on 100 Monte Carlo experiments.The best MSE value among different choices of γ is highlighted in bold.The bottom row shows estimated γ and the corresponding MSE and CI when estimated with our proposed method.The tuning parameters were set to (N, T, γ 0 ) = (2, 000, 300, 0.1) with 50 MCMC steps Experiment 3: Comparison withJewson and Rossell (2021)
Jewson and Rossell (2021)on and Rossell (2021), so we only provide an overview of the proof.Assume that, for simplicity, d y = 1.First one can show that H n (γ) can be decomposed into the sum of conditional expectation terms of H(y i , p γ (y i | θ, y −i )) = H(y i , p γ (y i | θ)), and the sum of conditional variance terms of