Linear Bayesian DRO

We discuss the DRO models based on misspecified parametric distribution assumptions. They are often imposed in newsvendor or portfolio loss with normal or exponential family; regression with gaussian assumptions.

Frequentist

Fit a single parametric nominal distribution to the observed data, then solve an f-DRO problem around that fitted distribution. Here we estimate a joint Gaussian distribution for (X, y) using its sample mean and maximum-likelihood covariance. We draw 500 joint samples from this fitted Gaussian and pass them to Chi2DRO.fit() for robust linear regression.

This is a finite Monte Carlo approximation to the parametric f-DRO problem. Increasing the simulated sample count improves the approximation; eps separately controls the chi-square ambiguity radius. The Gaussian parameters are point estimates, with no posterior sampling.

import numpy as np
from dro.linear_model.chi2_dro import Chi2DRO

# Observed regression data.
rng = np.random.default_rng(42)
X_observed = rng.normal(size=(100, 2))
y_observed = X_observed @ np.array([1.5, -2.0]) + rng.normal(scale=0.5, size=100)

# Fit one joint Gaussian nominal distribution and approximate it by sampling.
joint_observed = np.column_stack((X_observed, y_observed))
mean_hat = joint_observed.mean(axis=0)
cov_hat = np.cov(joint_observed, rowvar=False, ddof=0)
joint_nominal = rng.multivariate_normal(mean_hat, cov_hat, size=500)

frequentist_model = Chi2DRO(input_dim=2, model_type='ols', solver='CLARABEL')
frequentist_model.update({'eps': 0.1})
frequentist_params = frequentist_model.fit(joint_nominal[:, :-1], joint_nominal[:, -1])

print('Coefficients:', np.round(frequentist_params['theta'], 3))
print(f'Robust objective: {frequentist_model.robust_obj:.3f}')
print(f'Observed-data MSE: {frequentist_model.score(X_observed, y_observed):.3f}')

Coefficients: [ 1.443 -2.019]
Robust objective: 0.367
Observed-data MSE: 0.256

Bayesian

For the Bayesian-DRO model, we provide additional model type called newsvendor loss. The default setup is as follows: \(\ell(\theta;(X, Y)) = 3 (Y - \theta)_+ + (\theta - Y)_+\). In this case, X is the unit vector. That corresponds to exponential gamma conjugate.

That is, when conducting fitting, resampling, we ignore the existence of \(X\).

Both Bayesian examples explicitly set eps=0.1 to illustrate a positive-radius ambiguity set. The newsvendor example uses SCS for its exponential-cone formulation and fit_intercept=False, since its loss depends only on the order quantity.

## Newsvendor Example
import numpy as np
from dro.linear_model.chi2_dro import *
from dro.linear_model.bayesian_dro import *

from dro.data.dataloader_regression import regression_basic

feature_dim = 1
sample_num = 1000
np.random.seed(42)
X = np.ones((sample_num, feature_dim))
y = np.random.exponential(2, sample_num)

reg_model = BayesianDRO(
    input_dim=feature_dim, model_type='newsvendor',
    eps=0.1, solver='SCS', fit_intercept=False,
)
reg_model.update({'posterior_param_num': 1, 'distribution_class': 'Exponential', 'distance_type': 'KL'})
reg_model.fit(X, y)

{'theta': [3.5978596108504233], 'b': 0}

reg_model.distance_type

'KL'

## Regression Example
from dro.linear_model.chi2_dro import *
from dro.linear_model.bayesian_dro import *

from dro.data.dataloader_regression import regression_basic

feature_dim = 5
X, y = regression_basic(num_samples = 100, d = feature_dim, noise = 1)

reg_model = BayesianDRO(input_dim=feature_dim, model_type='ols', eps=0.1, solver='CLARABEL')

reg_model.update({'posterior_param_num': 5, 'distance_type': 'chi2'})
reg_model.fit(X, y)

sample lool 0
sample lool 1
sample lool 2
sample lool 3
sample lool 4

{'theta': [61.057566360599026,
  98.39982638076957,
  62.38125553686774,
  55.18121760481207,
  35.213916290192095],
 'b': array(-0.77980186)}