Distributionally Robust Optimisation with Bayesian Ambiguity Sets
Charita Dellaporta Patrick O'Hara Theodoros Damoulas

Bias correction of posterior means using MCMC outputs
    Yukito Iba
    September 9, 2024
=====================================================

*These authors contributed equally to this work.footnote
§ ABSTRACT 
 
Decision making under uncertainty is challenging since the data-generating process (DGP) is often unknown.
Bayesian inference proceeds by estimating the DGP through posterior beliefs about the model's parameters.
However, minimising the expected risk under these posterior beliefs can lead to sub-optimal decisions due to model uncertainty or limited, noisy observations.
To address this, we introduce Distributionally Robust Optimisation with Bayesian Ambiguity Sets (DRO-BAS) which hedges against uncertainty in the model by optimising the worst-case risk over a posterior-informed ambiguity set. We show that our method admits a closed-form dual representation for many exponential family members and showcase its improved out-of-sample robustness against existing Bayesian DRO methodology in the Newsvendor problem.





§ INTRODUCTION
 

Decision-makers are regularly confronted with the problem of optimising an objective under uncertainty. 
Let x ∈^d be a decision-making variable that minimises a stochastic objective function f: ^d ×Ξ→, where Ξ is the data space and let ^⋆∈(Ξ) be the data-generating process (DGP) where (Ξ) is the space of Borel distributions over Ξ. In practice, we do not have access to ^⋆ but to n independently and identically distributed (i.i.d.) observations := ξ_1:n∼^⋆. Without knowledge of the DGP, model-based inference considers a family of models _Θ := {_θ : θ∈Θ}⊂(Ξ) where each _θ has probability density function p(ξ | θ) for parameter space Θ⊆^k. In a Bayesian framework, data  is combined with a prior π(θ) to obtain posterior beliefs about θ through  Π(θ | ). Bayesian Risk Optimisation <cit.> then solves a stochastic optimisation problem:


    min_x ∈^d _θ∼Π(θ|)[_ξ∼_θ[f(x, ξ)]].


However, our Bayesian estimator is likely different from the true DGP due to model and data uncertainty: the number of observations may be small; the data noisy; or the prior or model may be misspecified.

The optimisation problem (<ref>) inherits any estimation error, and leads to overly optimistic decisions on out-of-sample scenarios even if the estimator is unbiased: this phenomenon is called the optimiser's curse <cit.>.
For example, if the number of observations is small and the prior is overly concentrated, then the decision is likely to be overly optimistic.

To hedge against the uncertainty of the estimated distribution, the field of Distributionally Robust Optimisation (DRO) minimises the expected objective function under the worst-case distribution that lies in an ambiguity set U ⊂(Ξ).
Discrepancy-based ambiguity sets contain distributions 
that are close to a nominal distribution in the sense of some discrepancy measure such as the Kullback-Leibler (KL) divergence <cit.>, Wasserstein distance <cit.> or Maximum Mean Discrepancy <cit.>.

For example, some model-based methods <cit.> consider a family of parametric models and create discrepancy-based ambiguity sets centered on the fitted model. 

However, uncertainty about the parameters is not captured in these
works, which can lead to a nominal distribution far away from the DGP when the data is limited. The established framework for capturing such uncertainty is Bayesian inference. 

The closest work to ours, using parametric Bayesian inference to inform the optimisation problem, is Bayesian DRO (BDRO) by <cit.>.

BDRO constructs discrepancy-based ambiguity sets with the KL divergence and takes an expected worst-case approach, under the posterior distribution.

More specifically, let U_ϵ(_θ) := {∈(Ξ) : (‖_θ) ≤ϵ} be the ambiguity set centered on distribution _θ with parameter ϵ∈ [0, ∞) controlling the size of the ambiguity set.

Under the expected value of the posterior, Bayesian DRO solves:


    min_x ∈^d _θ∼Π(θ|) [sup_∈ U_ϵ (_θ) _ξ∼[f(x, ξ)]],


where _θ∼Π(θ|)[Y] := ∫_Θ Y(θ)Π(θ|)  dθ denotes the expectation of random variable Y: Θ→ with respect to Π(θ|).

A decision maker is often interested in protecting against and quantifying the worst-case risk, but BDRO does not correspond to a worst-case risk analysis.
Moreover, the BDRO dual problem is a two-stage stochastic problem that involves a double expectation over the posterior and likelihood.
To get a good approximation of the dual problem, a large number of samples are required, which increases the solve time of the dual problem. 

We introduce DRO with Bayesian Ambiguity Sets (DRO-BAS), an alternative optimisation objective for Bayesian decision-making under uncertainty, based on a posterior-informed ambiguity set. The resulting problem corresponds to a worst-case risk minimisation over distributions with small expected deviation from the candidate model. We go beyond ball-based ambiguity sets, which are dependent on a single nominal distribution, by allowing the shape of the ambiguity set to be informed by the posterior.

For many exponential family models, we show that the dual formulation of DRO-BAS is an efficient single-stage stochastic program.



§ DRO WITH BAYESIAN AMBIGUITY SETS




We propose the following DRO-BAS objective:


    min_x ∈^d sup_: _θ∼Π[D(, _θ)] ≤ϵ _ξ∼ [ f_x(ξ) ],


where ∈(Ξ) is a distribution in the ambiguity set, f_x(ξ) := f(x, ξ) is the objective function, D: (Ξ) ×(Ξ) → is a divergence, and ϵ∈ [0, ∞) is a tolerance level.

The ambiguity set is informed by the posterior distribution Π by considering all probability measures ∈(Ξ) which are ϵ-away from _θ in expectation, with ϵ dictating the desired amount of risk in the decision. 

The shape of our ambiguity set is flexible and driven by the posterior distribution. This is contrary to standard ambiguity sets which correspond to a ball around a fixed nominal distribution. The DRO-BAS problem (<ref>) is still a worst-case approach, keeping with DRO tradition, instead of BDRO's expected worst-case formulation (<ref>), see <Ref>. 



The Bayesian posterior Π(θ|) targets the KL minimiser between the model family and ^⋆ <cit.>, hence it is natural to choose D(, _θ) to be the KL divergence of  with respect to _θ denoted by (‖_θ). This means that as n →∞ the posterior collapses to θ_0 := _θ∈Θ(^⋆, _θ) and the ambiguity set is just a KL-ball around _θ_0.


Using the KL divergence in the DRO-BAS problem in (<ref>), it is straight-forward to obtain an upper bound of the worst-case risk for general models (see <Ref> for a proof):

    
    sup_: _θ∼Π [(Q ‖_θ)] ≤ϵ _ξ∼[f_x(ξ)] ≤inf_γ≥ 0 γϵ + _θ∼Π[ γln_ξ∼_θ [ exp( f_x(ξ)/γ) ] ].



Exact closed-form solutions of DRO-BAS can be obtained for a wide range of exponential family models with conjugate priors. When the likelihood distribution is a member of the exponential family, a conjugate prior also belongs to the exponential family <cit.>.



In this setting, before we prove the main result, we start with an important Lemma.


 
    Let p(ξ|θ) be an exponential family likelihood and π(θ), Π(θ|) a conjugate prior-posterior pair, also members of the exponential family.
    Let τ_0, τ_n ∈ T be hyperparameters of the prior and posterior respectively, where T is the hyperparameter space.
    Let θ̅_n ∈Θ depend upon τ_n and let G: T → be a function of the hyperparameters. If the following identity holds:
    
    
    _θ∼Π[ ln p(ξ|θ) ] = ln p(ξ|θ̅_n) - G(τ_n),

    
    then the expected KL-divergence can be written as:
    
    
    _θ∼Π[  (‖_θ) ] =   (, _θ̅_n) + G(τ_n).



The condition in (<ref>) is a natural property of many exponential family models, some of which are showcased in Table <ref>. Future work aims to prove this for all exponential family models.

It is straightforward to establish the minimum tolerance level ϵ_min required to obtain a non-empty ambiguity set. Since the KL divergence is non-negative, under the condition of Lemma <ref>, for any ∈(Ξ):


    _θ∼Π [(‖_θ)] = (, _θ̅_n) + G(τ_n) ≥ G(τ_n) := ϵ_min(n).




We are now ready to prove our main result.
 
    Suppose the conditions of Lemma <ref> hold and ϵ≥ϵ_min(n) as in (<ref>).
    Let τ_n ∈ T, θ̅_n ∈Θ, and G: T →.
    Then
    
    sup_: _θ∼Π [(‖_θ)] ≤ϵ _ξ∼[f_x(ξ)] = inf_γ≥ 0 γ (ϵ - G(τ_n)) +  γln_ξ∼ p(ξ|θ̅_n) [ exp( f_x(ξ)/λ) ].



To guarantee that the DRO-BAS objective upper bounds the expected risk under the DGP, the decision-maker aims to choose ϵ large enough so that ^⋆ is contained in the ambiguity set. The condition in (<ref>) yields a closed-form expression for the optimal radius ϵ^⋆ by noting that:


    ϵ^⋆ = _θ∼Π [(^⋆‖_θ)] = (^⋆, _θ̅_n) + G(τ_n).


If the model is well-specified, and hence ^⋆ and _θ̅_n belong to the same exponential family, it is straightforward to obtain ϵ^⋆ based on the prior, posterior and true parameter values. 
We give examples in <Ref>.
In practice, since the true parameter values are unknown, we can approximate ϵ^⋆ using the observed samples. It follows that for any ϵ≥ϵ^⋆≥ϵ_min(n):


    _ξ∼^⋆[f(x, ξ)] ≤sup_: _θ∼Π [(Q ‖_θ)] ≤ϵ _ξ∼[f_x(ξ)].




§ THE NEWSVENDOR PROBLEM
 



  
Experiment setup. We evaluate DRO-BAS against the BDRO framework on a univariate Newsvendor problem with a well-specified univariate Gaussian likelihood with unknown mean and variance (<Ref> showcases a misspecified setting).
The goal is to choose an inventory level 0 ≤ x ≤ 50 of a perishable product with unknown customer demand ξ∈ that minimises the cost function f(x, ξ) = h max(0, x - ξ) + b max(0, ξ - x), where h and b are the holding cost and backorder cost per unit of the product respectively.
We let ^⋆ be a univariate Gaussian (μ_⋆, σ^2_⋆) with μ_⋆ = 25 and σ^2_⋆ = 100.
For random seed j = 1,…,200, the training dataset _n^(j) contains n = 20 observations and the test dataset _m^(j) contains m = 50 observations.
The conjugate prior and posterior are normal-gamma distributions (<Ref>).

N is the total number of samples from each model.
For each seed j, we run DRO-BAS and BDRO with N = 25, 100, 900 and across 21 different values of ϵ ranging from 0.05 to 3. For DRO-BAS, N is the number of samples from p(ξ|θ̅_n) and for BDRO, N = N_θ× N_ξ where N_θ is the number of posterior samples and N_ξ likelihood samples due to the double expectation present; we set N_θ = N_ξ to compare models on an equal N total samples regime. For a given ϵ, we calculate the out-of-sample mean m(ϵ) and variance v(ϵ) of the objective function f(x^(j)_ϵ, ξ̂_i) over all ξ̂_i ∈^(j)_m and over all seeds j=1,…,200, where x^(j)_ϵ is the optimal solution on training dataset ^(j)_n (see <Ref>).





  
Analysis. <Ref> shows that, for small sample size N = 25, 100, our framework dominates BDRO in the sense that DRO-BAS forms a Pareto front for the out-of-sample mean-variance tradeoff of the objective function f.
That is, for any ϵ_1, let m_BDRO(ϵ_1) and v_BDRO(ϵ_1) be the out-of-sample mean and variance respectively of BDRO:
then there exists ϵ_2 with out-of-sample mean m_BAS(ϵ_2) and variance v_BAS(ϵ_2) of BAS-DRO such that m_BAS(ϵ_2) < m_BDRO(ϵ_1) and v_BAS(ϵ_2) < v_BDRO(ϵ_1).
When N=900, <Ref> shows DRO-BAS and BDRO lie roughly on the same Pareto front.
To summarise, BDRO requires more samples N than DRO-BAS for good out-of-sample performance, likely because BDRO must evaluate a double expectation over the posterior and likelihood, whilst DRO-BAS only samples from p(ξ|θ̅_n).
For fixed N, the solve times for DRO-BAS and BDRO are broadly comparable (see <Ref>).







§ DISCUSSION


We proposed a novel approach to Bayesian decision-making under uncertainty through a DRO objective based on posterior-informed Bayesian ambiguity sets. The resulting optimisation problem is a single-stage stochastic program with closed-form formulation for a variety of exponential-family models. The suggested methodology has good out-of-sample performance, as showcased in <Ref>. Future work aims to extend DRO-BAS to a general formulation for exponential family models, including higher-dimensional problems, in which we expect to see further advantages of our method due to the nature of the Bayesian Ambiguity Set. 


CD acknowledges support from EPSRC grant [EP/T51794X/1] as part of the Warwick CDT in Mathematics and Statistics. PO and TD acknowledge support from a UKRI Turing AI acceleration Fellowship [EP/V02678X/1] and a Turing Impact Award from the Alan Turing Institute. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) license to any Author Accepted Manuscript version arising from this submission.





  

    

    Supplementary Material
    

  



The Supplementary Material is organised as follows:  <Ref> provides details for all the exponential family models discussed in <Ref>, while <Ref> contains the proofs of all mathematical results appearing in the main text. <Ref> provides additional experimental details for the Newsvendor problem in <Ref>. Finally, in <Ref> we present experimental results for the newsvendor problem example with a misspecified model. 



§ SPECIAL CASES
 

We derive the values of G(τ_n) and θ̅_n for different likelihoods and conjugate prior/posterior in <Ref>.
Each subsection contains a corollary with the result in <Ref>.




 §.§ Gaussian model with unknown mean and known variance


Let the random variable ξ be univariate and have continuous support.
We assume the variance σ^2 of ξ is known.
We estimate the mean of a univariate Gaussian distribution with known variance σ^2.
The example can be found in <cit.>.
We define our parameter θ to be the unknown mean μ and we place a Gaussian prior π(μ) over it.

The likelihood is p(ξ|μ) = (μ, σ^2), the prior over μ is 
π(μ) = (μ_0, σ_0^2), and the conjugate posterior is Π(μ|) = (μ|μ_n , σ_n^2), 
where

    μ_n := σ^2/nσ^2_0 + σ^2μ_0 + n σ_0^2/n σ_0^2 + σ^2μ̂      1/σ^2_n := 1/σ_0^2 + n/σ^2      μ̂ := 1/n∑_i=1^n ξ_i.



 
    Let μ_n ∈ and σ, σ_n ∈_+. Then
    
    _μ∼(μ_n,σ_n^2)[ log(μ, σ^2) ] = log(μ_n, σ^2) - σ_n^2/2σ^2.



    The result is a special case of <cit.>.


 
    When the likelihood is a Gaussian distribution with unknown mean and known variance σ^2 and the prior and posterior are Gaussian distributions (see <Ref>), then <Ref> holds with θ̅_n = μ_n and G(τ_n) = σ_n^2/2σ^2.


        <Ref> shows that the condition (<ref>) in <Ref> holds, thus <Ref> follows.




  
Tolerance level ϵ 

In the well-specified case, where we assume that ^⋆ := _θ^⋆ for some θ^⋆∈Θ, it is easy to obtain the required size of the ambiguity set exactly. Let θ^⋆ := μ^⋆ and ^⋆ := N(μ^⋆, σ^2). By <Ref> it follows that:


    _θ∼ N(μ_n, σ^2_n)[(^⋆, _θ) ] = (μ^⋆ - μ_n)^2 + σ^2_n/2 σ^2.


So for a fixed finite sample ξ_1:n∼^⋆, if ϵ≥ϵ^⋆ := (μ^⋆ - μ_n)^2 + σ^2_n/2 σ^2, it follows that DRO-BAS upper-bounds the target optimisation objective: 


    _ξ∼^⋆ [f_x(ξ)] ≤sup_: _θ∼Π(·|ξ_1:n) [(Q ‖ _θ)] ≤ϵ _ξ∼[f_x(ξ)].
 

In practice, since μ^⋆ is unknown, ϵ^⋆ can be approximated by using the sample mean. 



 §.§ Gaussian Model with Unknown Mean and Variance
 

In this section, we consider a Bayesian model that estimates the unknown mean and variance of a uni-variate Gaussian distribution.
We define our model in <Ref>, then prove some preliminary results for the normal-gamma distribution, before proving our main result in <Ref>.


Following <cit.>, we place a normal-gamma prior over the mean μ and precision λ = σ^-2.
The normal-gamma prior is the conjugate to a Gaussian likelihood and results in a normal-gamma posterior distribution. 

    Likelihood:       p(ξ|μ, λ) = (ξ|μ, λ^-1) 
       Prior:      π(μ, λ) = NG(μ, λ|μ_0, κ_0, α_0, β_0 ) = (μ|μ_0, (κ_0 λ)^-1) · Ga(λ|α_0, β_0) 
       Posterior:      π(μ, λ|) = NG(μ, λ|μ_n, κ_n, α_n, β_n) = (μ|μ_n , (λκ_n)^-1 ) · Ga(λ|α_n, β_n)

where

    μ_n := κ_0 μ_0 + n ξ̅_n/n + κ_0, κ_n := κ_0 + n, α_n := α_0 + n/2, ξ̅_n = 1/n∑_i-1^n ξ_i,
       β_n := β_0 + 1/2∑_i=1^n (ξ_i - ξ̅_n)^2 + κ_0 n (ξ̅_n - μ_0)^2/2(κ_0 + n).



In <Ref>, we derive condition <Ref> for a Gaussian model with unknown mean and unknown variance.
Before proceeding, we need to define the gamma and digamma functions and recall the moments of the normal-gamma distribution.

 
    The gamma function Γ: → and digamma function ψ: → are
    
    Γ(z) := (z-1)!       ψ(z) := /ẓlnΓ(z).



 
    Let NG(μ, λ|μ_n, κ_n, α_n, β_n) be a normal-gamma distribution with parameters μ_n ∈ and κ_n, α_n, β_n ∈_+.
    The moments of the normal-gamma distribution are
    
    _NG[lnλ] = ψ(α_n) - lnβ_n,       _NG[λ] = α_n/β_n,       _NG[λμ] = μ_n α_n/β_n,       _NG[λμ^2] = 1/κ_n + μ_n^2α_n/β_n.

    where ψ: → is the digamma function from <Ref>.


 
    Let μ_n ∈ and κ_n, α_n, β_n ∈_+. Then
    
    _NG(μ, λ|μ_n, κ_n, α_n, β_n)[ ln(ξ|μ, λ^-1 ) ] = ln( ξ|μ_n, β_n/α_n) - 1/2( 1/κ_n +lnα_n - ψ(α_n) )

    where ψ: → is the digamma function from <Ref>.


    First, observe that the natural logarithm of the Gaussian distribution may be re-written as
    
    ln(ξ|μ, λ^-1)    = ln( λ^1/2/(2π)^1/2exp( - λ/2 (ξ - μ)^2 ) ) 
       = 1/2lnλ - 1/2ln 2π - λ/2 (ξ - μ)^2 
       = 1/2lnλ - 1/2ln 2π - 1/2λξ^2 + λμξ - 1/2λμ^2.

    In what follows, for shorthand, we denote the expectation _NG(μ, λ|μ_n, κ_n, α_n, β_n) as _μ, λ∼ NG:
    
    
    _μ, λ∼ NG[ ln(ξ|μ, λ^-1 ) ] 
         (i)=_μ, λ∼ NG[ 1/2lnλ - 1/2ln 2π - 1/2λξ^2 + λμξ - 1/2λμ^2 ] 
         (ii)= - 1/2ln 2π + 1/2_μ, λ∼ NG[ lnλ] - 1/2ξ^2 ·_μ, λ∼ NG[ λ] + ξ·_μ, λ∼ NG[ λμ] - 1/2_μ, λ∼ NG[ λμ^2 ] 
         (iii)= - 1/2ln 2π + 1/2 (ψ(α_n) - lnβ_n) - 1/2ξ^2 α_n/β_n + ξμ_n α_n/β_n - 1/2(1/κ_n + μ_n^2α_n/β_n) 
         (iv)= - 1/2ln 2π + 1/2(ψ(α_n) - lnβ_n - 1/κ_n) - α_n/2β_n( ξ - μ_n )^2 
         (v)= - 1/2ln 2π - 1/2(lnβ_n - lnα_n + lnα_n - ψ(α_n) + 1/κ_n) + lnexp(- 1/2β_n/α_n( ξ - μ_n )^2 )
         (vi)= - 1/2(lnα_n - ψ(α_n) + 1/κ_n) - 1/2ln(2πβ_n/α_n) + lnexp(- 1/2β_n/α_n( ξ - μ_n )^2 ) 
         (vii)= - 1/2( 1/2α_n + I_α_n + 1/κ_n) + ln( 1/√(2πβ_n/α_n)exp( - 1/2 β_n/α_n(ξ - μ_n)^2 ) ) 
         (viii)= - 1/2( lnα_n - ψ(α_n) + 1/κ_n) + ln( ξ|μ_n, β_n/α_n)

    where in equation (i) we take the expectation over the normal-gamma distribution;
    (ii) we apply linearity of expectation;
    (iii) we use the moment-generating functions from <Ref>;
    (iv) we complete the square;
    (v) we add and subtract lnα_n;
    (vi) and (vii) we re-arrange and apply log identities;
    and finally in (viii) we use the definition of a Gaussian probability density function.


 
    When the likelihood is a Gaussian distribution with unknown mean and variance, and the conjugate prior and posterior are normal-gamma distributions (see <Ref>), then <Ref> holds with θ̅_n = (μ_n, β_n/α_n) and G(τ_n) = 1/2( lnα_n - ψ(α_n) + 1/κ_n).


    <Ref> shows that the condition (<ref>) in <Ref> holds, thus <Ref> follows.




  
Tolerance level ϵ

In the well-specified case, where we assume that ^⋆ := _θ^⋆ for some θ^⋆∈Θ, it is easy to obtain the required size of the ambiguity set exactly. Let θ^⋆ := (μ^⋆, λ^⋆^-1) and ^⋆ := N(μ^⋆, λ^⋆^-1). Using <Ref> we obtain:


    _μ, λ∼ NG(μ, λ|μ_n, κ_n, α_n, β_n)[(^⋆, (ξ|μ, λ^-1)) ] 
            = ( ^⋆ ‖ (μ_n, β_n/α_n) ) + 1/2( 1/κ_n + lnα_n - ψ(α_n) ) 
            = ln(√(λ^⋆β_n/α_n)) + λ^⋆^-1 + (μ^⋆ - μ_n)^2/2 β_n/α_n - 1/2 + 1/2( 1/κ_n + lnα_n - ψ(α_n) ) 
            = 1/2(ln(λ^⋆β_n )
        + λ^⋆^-1 + (μ^⋆ - μ_n)^2)/2 β_n/α_n - 1
        + 1/κ_n - ψ(α_n)
        ).





 §.§ Exponential likelihood with conjugate gamma prior


 
   The likelihood p(ξ|θ) is an exponential distribution (ξ|θ) where θ > 0 is the rate parameter.
    The prior π(θ) is a gamma distribution (θ|α_0, β_0) with shape α_0 > 0 and rate β_0 > 0.
    The parameters α_n, β_n of the posterior π(θ|) = (θ|α_n, β_n) are given by α_n = α_0 + n and β_n = β_0 + ∑_ξ_i ∈ ξ_i. 



 
    When the likelihood is an exponential distribution with gamma prior and posterior (see <Ref>), then
    
    _(θ|α_n, β_n)[ ln(ξ|θ) ] = ln( ξ|α_n/β_n) + ψ(α_n) - lnα_n.



    Starting from the left-hand side, we take the log of the PDF of the exponential distribution, then use the logarithm expectation of the gamma distribution, and finally re-arrange using log identities:
    
    _(θ|α_n, β_n)[ ln(ξ|θ) ]    = _(θ|α_n, β_n)[ lnθ - θξ] 
       = ψ(α_n) - lnβ_n - α_n/β_nξ
       = ψ(α_n) - lnα_n + lnα_n/β_n - α_n/β_nξ
       = ψ(α_n) - lnα_n + ln( α_n/β_nexp( - α_n/β_nξ) ) 
       = ψ(α_n) - lnα_n + ln(ξ|α_n/β_n).

    The last line follows by the definition of the PDF of the exponential distribution.


 
    When the likelihood is an exponential distribution with gamma prior and posterior, then <Ref> holds with θ̅_n = α_n/β_n and G(τ_n) = ψ(α_n) - lnα_n.


    <Ref> shows that the condition (<ref>) in <Ref> holds, thus <Ref> follows.




  
Tolerance level ϵ

In the well-specified case, where we assume that ^⋆ := _θ^⋆ for some θ^⋆∈Θ, it is easy to obtain the required size of the ambiguity set exactly. Let θ^⋆ be the true rate parameter, i.e. ^⋆ := (θ^⋆). Using <Ref> we obtain:


    _(θ|α_n, β_n)[(^⋆, (θ) ] 
            = ( ^⋆ ‖ (α_n/β_n) ) + ψ(α_n) - ln(α_n)  
            = ln(θ^⋆) - ln(α_n/β_n) + α_n/β_n θ^⋆ - 1 + ψ(α_n) - ln(α_n).




§ PROOFS OF THEORETICAL RESULTS
 


 §.§ Proofs of DRO-BAS upper bound in Equation (<ref>)
 

Before proving the required upper bound, we recall the definition of the KL divergence and its convex conjugate. 


    Let μ,ν∈(Ξ) and assume μ is absolutely continuous with respect to ν (μ≪ν).
    The KL-divergence of μ with respect to ν is defined as: 
    
    (μ‖ν):=  ∫_Ξln(μ(dξ)/ν(dξ)) μ(dξ).




    Let ν∈(Ξ) be non-negative and finite.
    The convex conjugate ^⋆(·‖ν) of (·‖ν) is
    
    ^⋆(·‖ν)(h) = ln( ∫_Ξexp(h) dν).



    See Proposition 28 and Example 7 in <cit.>.





    The result follows from a standard Lagrangian duality argument and an application of Jensen's inequality.
    More specifically, we introduce a Lagrangian variable γ≥ 0 for the expected-ball constraint on the left-hand side of (<ref>) as follows:
    
    sup_: _θ∼Π [(Q ‖_θ)] ≤ϵ _Q[f_x]
               (i)≤inf_γ≥ 0 sup_∈(Ξ) _Q[f_x] + γϵ - γ_Π[  (‖_θ) ] 
       (ii)=inf_γ≥ 0 γϵ + sup_∈(Ξ) _Q[f_x] -  _Π[ γ(‖_θ)  ] 
       (iii)=inf_γ≥ 0 γϵ + ( _Π[  γ(·‖_θ) ] )^⋆ (f_x) 
       (iv)≤inf_γ≥ 0 γϵ + _Π[ ( γ(·‖_θ) )^⋆(f_x)  ] 
       (v)=inf_γ≥ 0 γϵ + _Π[ γln__θ[ exp(f_x/γ) ] ].

    Inequality (i) holds by weak duality.
    Equality (ii) holds by linearity of expectation and a simple rearrangement.
    Equality (iii) holds by the definition of the conjugate function.
    Inequality (iv) holds by Jensen's inequality ([·])^⋆≤[(·)^⋆] because the conjugate is a convex function.
    Equality (v) holds by <Ref> and the fact that for γ≥ 0 and function ϕ, (γϕ)^⋆(y) = γϕ^⋆(y/γ).




 §.§ Proof of <Ref>

Starting from the left-hand side, we have
    
    _Π[  (‖_θ) ]
           (i)=_θ∼π(θ|)[ ∫_Ξ q(ξ) ln( q(ξ)/p(ξ|θ)) ξ̣] 
       (ii)=_θ∼π(θ|)[ ∫_Ξ q(ξ) ln( q(ξ) ) - q(ξ) ln( p(ξ|θ)) ξ̣] 
       (iii)=∫_Ξ q(ξ) ln( q(ξ) ) - q(ξ) ·_θ∼π(θ|)[ ln( p(ξ|θ)) ] ξ̣
       (iv)=∫_Ξ q(ξ) ln( q(ξ) ) - q(ξ) ·(ln p(ξ|θ̅_n) - G(τ_n) ) ξ̣
       (v)=∫_Ξ q(ξ) ln( q(ξ)/p(ξ|θ̅_n)) ξ̣ + ∫_Ξ  q(ξ) · G(τ_n) ξ̣
       (vi)=(‖_θ̅_n) + _[G(τ_n)] 
       (vii)=(q(ξ) ‖ p(ξ|θ̅_n)) + G(τ_n).

where (i) is by the definition of the KL-divergence; (ii) follows by log properties; (iii) holds by linearity of expectation; (iv) holds by condition (<ref>) in <Ref>; (v) holds by rearrangement and properties of log; (vi) holds by the definition of the KL-divergence and the definition of _; and (vii) holds by the expected value of a constant.




 §.§ Proof of <Ref>



We begin by restating the Lagrangian dual from the proof of <Ref>, but with the added claim that strong duality holds between the primal and dual problems:

    sup_: _θ∼Π [(‖_θ)] ≤ϵ _[f_x] =inf_γ≥ 0 sup_∈(Ξ) _[f_x] + γϵ - γ_Π[  (‖_θ) ].




The conditions under which our claim of strong duality holds will be proved later.
Next, we substitute the right-hand side of equation (vi) above into the dual problem in (<ref>):

    sup_: _θ∼Π [( ‖ _θ)] ≤ϵ _ξ∼[f_x(ξ)] 
            = inf_γ≥ 0 γϵ + sup_∈(Ξ) ∫_Ξ f_x(ξ) q(ξ) ξ̣ - γ( ( ‖  p(ξ|θ̅_n) ) + G(τ_n) ) 
           = inf_γ≥ 0 γϵ + γ G(τ_n) + ( γ (· ‖  p(ξ|θ̅_n) ) )^⋆( f_x(ξ) ) 
           = inf_γ≥ 0 γ (ϵ - G(τ_n)) + γln_ξ∼  p(ξ|θ̅_n)[ exp( f_x(ξ)/γ) ],

where the second and third equality holds by the definition of the conjugate of the KL-divergence and by <Ref>.

Finally, it remains to argue that strong duality holds. First, note that the primal problem is a concave optimisation problem with respect to distribution . Second, when ϵ > G(τ_n), then distribution = p(ξ|θ̅_n) is a strictly feasible point to the primal constraint because

    _θ∼Π [( ‖ _θ)] = 0 < ϵ - G(τ_n).





§ NEWSVENDOR PROBLEM - ADDITIONAL DETAILS


We provide additional details about our Newsvendor experiment in <Ref> when ^⋆ is a Gaussian distribution with μ_⋆ = 25 and σ^2_⋆ = 100.



  
Hyperparameters. The prior and posterior are normal-gamma distributions.
We set the prior hyperparameters to be μ_0 = 0 and κ_0, α_0, β_0 = 1.
The derivation of the hyperparameters can be found in <Ref>.



  
Values of ϵ_min and ϵ^⋆. From <Ref> and equation (<ref>), the value of ϵ_min is 0.047.
From equation (<ref>), the average value of ϵ^⋆ over all J seeds is 0.089 with standard deviation 0.048.



  
Implementation. We implemented the dual problems for DRO-BAS (<Ref>) and BDRO <cit.> in Python using CVXPY version 1.5.2 and the MOSEK solver version 10.1.28. Our implementation uses disciplined parametrized programming <cit.> which - after an initial warm start for seed j=1 - allows us to solve subsequent seeds j = 2,…,J rapidly (see <Ref>).
We used a 12-core Dual Intel Xeon E5-2643 v3 @ 3.4 Ghz with 128GB RAM.



  
Out-of-sample mean and variance. For a given ϵ and seed j, let the optimal solution be x^(j)(ϵ).
We calculate the out-of-sample mean m^(j)(ϵ) = _ξ∼^(j)_m [f(x^(j)(ϵ),ξ)] and variance v^(j)(ϵ) = Var_ξ∼^(j)_m[ f(x^(j)(ϵ),ξ) ] of the objective under the empirical test distribution ^(j)_m.
For a given ϵ, the out-of-sample mean m(ϵ) and variance v(ϵ) across all seeds is

    m(ϵ) = 1/J∑_j=1^J m^(j)(ϵ)        v(ϵ) = 1/J∑_j=1^J v^(j)(ϵ) + 1/J-1∑_j=1^J (m^(j)(ϵ) - m(ϵ) )^2.

The out-of-sample variance v(ϵ) is equal to the mean of the variances v^(j)(ϵ) plus the variance of the means m^(j)(ϵ) <cit.>.





  
Solve time. On the initial warm-start seed j=1, for each N, DRO-BAS solves the dual problem from <Ref> faster than the BDRO dual problem.
For example, when N=900, DRO-BAS solves problems in 0.27 seconds on average, whilst BDRO solves problem in 5.56 seconds.
These results suggest that, for fixed N, if the solve is started from scratch with no warm start, then DRO-BAS will solve instances faster than BDRO.
For seeds j=2,…,J, disciplined parametrized programming (DPP) significantly speeds up the solve for BDRO: the average solve time for seeds j=2,…,J is 0.40 seconds when N=900.
In contrast, DPP does not speed up the solve for DRO-BAS: the average solve time for seeds j=2,…,J is 0.40 when N=900.
We conjecture that the speed up for BDRO using DPP is because BDRO has N_θ Lagrangian dual variables compared to DRO-BAS having exactly one Lagrangian dual variable.
BDRO then benefits from the warm start because it can reuse the presolve effort spent on the N_θ dual variables spent during the warm start.






§ MISSPECIFIED SUPPLEMENTARY EXPERIMENTS - TRUNCATED NORMAL
 

In this section, we present additional experiments when the data-generating process ^⋆ is a truncated normal distribution.
The truncated normal has mean μ_⋆ = 10 and variance σ^2_⋆ = 100.
The likelihood is a Gaussian distribution, so our model is misspecified.
The conjugate prior and posterior are still normal-gamma distributions with the same hyperparameters as <Ref>.
The experimental setup is also the same as <Ref>: the values of ϵ, N, n, m, and J are all specified the same.



  
Analysis. When the likelihood is misspecified, <Ref> shows the out-of-sample mean-variance tradeoff is again a Pareto front.
This is the same conclusion as the well-specified case in <Ref>.
Furthermore, when N=900, DRO-BAS has a small advantage on the mean-variance tradeoff.



  
Solve time. <Ref> shows the same conclusions about the solve time from <Ref> can be made about the solve time for the truncated normal data-generating process.




