259  Robustness to Natural Distribution Shift

259.1 1. Natural Versus Adversarial Distribution Shift

Machine learning systems deployed in the real world encounter a fundamental challenge: the statistical distribution of inputs at inference time rarely matches the distribution seen during training. This mismatch is called distribution shift, and it manifests in two fundamentally distinct forms.

Adversarial shift is worst-case, imperceptible, and intentional. An attacker with knowledge of the model crafts a perturbation \(\delta\) with \(\|\delta\|_p \leq \epsilon\) that maximizes loss while remaining visually indistinguishable from the original input. The perturbation is tailored to exploit a specific model’s decision boundary and does not generalize across models or occur spontaneously in nature.

Natural distribution shift is benign, unintentional, and pervasive. It arises from real-world variation: changes in lighting, weather, camera optics, geographic location, demographic composition of users, or simply the passage of time. A model trained on summer street photographs may degrade on winter images not because an adversary crafted the snow, but because snowflakes were absent from the training set. A clinical model trained on data from one hospital system may perform poorly at another because staining protocols differ. These shifts are benign in intent but catastrophic in consequence when models are assumed to be robust by default.

Both forms cause accuracy degradation, but natural shift is far more prevalent in deployment. Adversarial robustness and natural robustness are often studied separately and have different solutions: adversarial training against \(\ell_p\) perturbations does not reliably improve natural corruption robustness, and augmentation-based methods that improve natural robustness do not confer adversarial guarantees (Ford et al. 2019). The present chapter focuses entirely on natural distribution shift.

The mathematical framing is as follows. Let \(P_\text{train}\) denote the training distribution over \((X, Y)\) pairs and let \(\{P_e\}_{e \in \mathcal{E}}\) denote a collection of environments or domains that arise at test time. Empirical risk minimization (ERM) finds parameters \(\theta\) minimizing

\[\hat{\theta}_\text{ERM} = \argmin_\theta \mathbb{E}_{(x,y) \sim P_\text{train}} [\ell(f_\theta(x), y)]\]

but offers no guarantee about performance under \(P_e\) for \(e \notin \text{support}(P_\text{train})\). The goal of robust machine learning is to produce models whose worst-case or average-case performance across the \(P_e\) is high, without access to labeled target-domain data at training time.

259.2 2. Benchmark Suites for Natural Robustness

Rigorous empirical measurement of natural robustness requires controlled benchmark suites. Several have become standard.

259.2.1 2.1 ImageNet-C

Hendrycks and Dietterich (2019) introduced ImageNet-C to enable systematic evaluation of robustness to common visual corruptions. The benchmark applies 15 algorithmically defined corruption types to the ImageNet validation set at 5 severity levels each, yielding 75 distinct corrupted datasets. The corruption types span four broad categories: noise (Gaussian, shot, impulse), blur (defocus, glass, motion, zoom), weather (snow, frost, fog, brightness), and digital artifacts (contrast, elastic transform, pixelate, JPEG compression).

Performance is summarized by the mean Corruption Error (mCE):

\[\text{mCE} = \frac{1}{15} \sum_{c=1}^{15} \frac{\sum_{s=1}^{5} \text{CE}_{c,s}^f}{\sum_{s=1}^{5} \text{CE}_{c,s}^{\text{AlexNet}}}\]

where \(\text{CE}_{c,s}^f\) is the top-1 error of model \(f\) on corruption \(c\) at severity \(s\), normalized by AlexNet’s error on the same corruption. An mCE below 100 indicates better-than-AlexNet robustness; state-of-the-art models trained with strong augmentation achieve mCE near 40.

ImageNet-C has revealed that accuracy on clean ImageNet does not predict corruption robustness. A model with 3% lower clean error than a baseline often shows worse mCE if it was not trained with robustness-promoting augmentation.

259.2.2 2.2 ImageNet-R and ImageNet-A

Hendrycks et al. (2021) introduced two complementary benchmarks targeting different aspects of distribution shift.

ImageNet-R (Rendition) contains 30,000 images spanning 200 ImageNet classes, collected from sources such as art, cartoons, graffiti, embroidery, origami, sculptures, and toys. The labels are identical to ImageNet, but the visual style is radically different. ImageNet-R tests whether models learn object shape and structure rather than texture and context: Geirhos et al. (2019) showed that standard CNNs are strongly texture-biased, which explains large accuracy drops when the rendering style changes.

ImageNet-A (Adversarial) contains 7,500 naturally occurring images that were selected specifically because a ResNet-50 misclassifies them, achieving 0% top-1 accuracy. These are not adversarially perturbed; they are real photographs (e.g., atypical viewpoints, cluttered backgrounds, unusual color distributions) that happen to be highly challenging for standard models. Accuracy on ImageNet-A correlates strongly with model size and robustness training, making it a useful stress test.

259.2.3 2.3 WILDS

Koh et al. (2021) introduced the WILDS benchmark to address a gap in prior work: ImageNet-C and ImageNet-R measure shifts under controlled corruption or style change, but real deployment shifts are harder to characterize. WILDS compiles 10 real-world datasets spanning biology, medicine, NLP, and vision, each with documented distribution shifts that occur in practice.

Examples include: iWildCam (wildlife cameras across geographic locations and time), Camelyon17 (histology slides from five hospitals), Amazon (product reviews across time and user demographics), and Poverty (satellite imagery across African countries). Each dataset includes source domains for training and target domains for evaluation, with no target-domain labels available during training. WILDS exposes that methods effective on synthetic shifts (ImageNet-C) sometimes fail on real-world shifts, motivating techniques that operate on genuine distributional heterogeneity.

259.3 3. Why Empirical Risk Minimization Fails Under Shift

Standard training via ERM is vulnerable to natural distribution shift for a well-understood reason: it minimizes average loss over the training distribution, which allows it to exploit any statistical regularity that reduces training loss, regardless of whether that regularity reflects the underlying causal structure.

Consider a model trained to classify cows versus camels. If the training set contains cows predominantly on grass backgrounds and camels predominantly on desert backgrounds, ERM can achieve near-zero training error by learning a texture-based rule: “grass-like background \(\Rightarrow\) cow.” This spurious correlation (Sagawa et al. 2020) is sufficient for training accuracy but breaks immediately when a cow appears on a beach or a camel appears in a zoo enclosure with grass.

Formally, decompose the input features into invariant features \(Z^I\) that have a stable causal relationship to \(Y\) across environments, and spurious features \(Z^S\) whose correlation with \(Y\) is environment-dependent. ERM cannot distinguish \(Z^I\) from \(Z^S\) using data from a single training distribution. It will incorporate \(Z^S\) into its decision rule whenever doing so reduces training loss.

This failure mode is exacerbated by the over-parameterization of modern deep networks: a sufficiently flexible model can memorize arbitrary dataset-specific patterns, including background textures, acquisition artifacts, and demographic correlates that are not causally linked to the label. The resulting model may have high accuracy on the training distribution and its i.i.d. test set but fail catastrophically under any shift that changes the distribution of \(Z^S\) while leaving \(Y\) unchanged.

The solution space branches into two directions: (1) data augmentation to make the training distribution broader, reducing the statistical advantage of spurious features; and (2) explicit regularization to encourage the learned representation to be invariant across environments. We treat each in turn.

259.4 4. Augmentation Strategies for Corruption Robustness

259.4.1 4.1 AugMix

Hendrycks et al. (2020) proposed AugMix, which improves corruption robustness without requiring corrupted training data. The algorithm stochastically samples chains of augmentation operations, applies them to the input image to obtain several augmented views \(x_\text{aug}^{(k)}\), and mixes them with the original:

\[x_\text{mixed} = \lambda x_\text{orig} + (1-\lambda) \sum_{k} w_k x_\text{aug}^{(k)}\]

where \(\lambda \sim \text{Beta}(\alpha, \alpha)\) and weights \(w_k\) are sampled from a Dirichlet distribution. The augmentation chains are drawn from a set including rotation, translation, shear, solarize, equalize, autocontrast, and color jitter, applied in random order and random magnitude.

The key innovation is a Jensen-Shannon consistency loss applied alongside the standard cross-entropy loss. For the original image and two independently augmented views \(x_1, x_2\), the model predicts distributions \(p_0, p_1, p_2\) and the additional loss term is

\[\mathcal{L}_\text{JS} = \frac{1}{3}\left[D_\text{KL}(p_0 \| M) + D_\text{KL}(p_1 \| M) + D_\text{KL}(p_2 \| M)\right], \quad M = \frac{p_0 + p_1 + p_2}{3}\]

This encourages the model to produce consistent predictions across different views, which has the effect of smoothing the decision boundary with respect to augmentation-style variations. AugMix reduces mCE on ImageNet-C by approximately 12 percentage points relative to standard training without accuracy loss on clean ImageNet.

259.4.2 4.2 RandAugment

Cubuk et al. (2020) proposed RandAugment as a simplified alternative to learned augmentation policies. At each training step, \(N\) augmentation operations are selected uniformly at random from a predefined library of 14 (including Translate, Rotate, AutoContrast, Equalize, Solarize, Color, Posterize, Contrast, Brightness, Sharpness, ShearX, ShearY, CutOut). All selected operations are applied at the same global magnitude \(M\), which is a single hyperparameter tunable by validation. RandAugment eliminates the need for a separate proxy task to learn per-operation magnitudes, making it practically simpler while achieving competitive performance.

259.4.3 4.3 DeepAugment

Hendrycks et al. (2021) introduced DeepAugment, which uses image-to-image neural networks (such as EDSR super-resolution and a U-Net) as augmentation engines. Inputs are passed through these networks with their weights randomly perturbed, producing extreme but semantically valid transformations that are difficult to specify analytically. Combined with AugMix, DeepAugment achieves state-of-the-art robustness on ImageNet-C among models without extra training data.

259.5 5. Test-Time Augmentation

Test-time augmentation (TTA) is a simple, zero-cost-at-train-time technique for improving robustness at inference: rather than forwarding a single test image, forward \(T\) augmented versions and average the softmax predictions:

\[\hat{p}(y|x) = \frac{1}{T} \sum_{t=1}^{T} f_\theta(g_t(x))\]

where \(g_t\) are random augmentation functions (crops, flips, color jitter). Averaging over multiple views reduces variance in the prediction and is particularly effective when the model’s predictions are inconsistent across augmentations, which is precisely what happens under distribution shift. TTA requires no retraining and scales trivially with compute.

import torch
import torch.nn as nn
import torchvision
import torchvision.transforms as T
from torch.utils.data import DataLoader

torch.manual_seed(0)

# Small CNN for CIFAR-10
class SmallCNN(nn.Module):
    def __init__(self):
        super().__init__()
        self.features = nn.Sequential(
            nn.Conv2d(3, 32, 3, padding=1), nn.ReLU(),
            nn.Conv2d(32, 64, 3, padding=1), nn.ReLU(),
            nn.MaxPool2d(2),
            nn.Conv2d(64, 128, 3, padding=1), nn.ReLU(),
            nn.MaxPool2d(2),
        )
        self.classifier = nn.Sequential(
            nn.Flatten(),
            nn.Linear(128 * 8 * 8, 256), nn.ReLU(),
            nn.Linear(256, 10),
        )

    def forward(self, x):
        return self.classifier(self.features(x))


# Minimal training loop (few epochs for demonstration)
mean = (0.4914, 0.4822, 0.4465)
std  = (0.2023, 0.1994, 0.2010)

train_tf = T.Compose([T.RandomHorizontalFlip(), T.ToTensor(), T.Normalize(mean, std)])
test_tf  = T.Compose([T.ToTensor(), T.Normalize(mean, std)])

train_ds = torchvision.datasets.CIFAR10(root="/tmp/cifar", train=True,  download=True, transform=train_tf)
test_ds  = torchvision.datasets.CIFAR10(root="/tmp/cifar", train=False, download=True, transform=test_tf)

train_loader = DataLoader(train_ds, batch_size=256, shuffle=True,  num_workers=0)
test_loader  = DataLoader(test_ds,  batch_size=256, shuffle=False, num_workers=0)

model = SmallCNN()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
criterion = nn.CrossEntropyLoss()

for epoch in range(3):
    model.train()
    for x, y in train_loader:
        optimizer.zero_grad()
        criterion(model(x), y).backward()
        optimizer.step()

# Standard evaluation
model.eval()
correct = total = 0
with torch.no_grad():
    for x, y in test_loader:
        preds = model(x).argmax(1)
        correct += (preds == y).sum().item()
        total += len(y)
print(f"Clean accuracy: {100 * correct / total:.1f}%")

# TTA evaluation: average over T random crop + flip views
def tta_predict(model, x, T=8):
    aug = T_tta = T.Compose([
        T.RandomHorizontalFlip(),
        T.RandomCrop(32, padding=4),
    ])
    preds = torch.zeros(x.shape[0], 10)
    for _ in range(T):
        x_aug = torch.stack([aug(img) for img in x])
        with torch.no_grad():
            preds += torch.softmax(model(x_aug), dim=1)
    return preds.argmax(1)

tta_correct = tta_total = 0
for x, y in test_loader:
    preds = tta_predict(model, x, T=8)
    tta_correct += (preds == y).sum().item()
    tta_total += len(y)
print(f"TTA accuracy (T=8): {100 * tta_correct / tta_total:.1f}%")

259.6 6. Measuring the Natural Robustness Gap

The robustness gap is the difference in accuracy between the clean test set and a corrupted test set. The following code trains a small CNN on clean CIFAR-10 and measures accuracy as a function of Gaussian noise corruption severity, demonstrating how rapidly standard models degrade.

import torch
import torch.nn as nn
import torchvision
import torchvision.transforms as T
from torch.utils.data import DataLoader, TensorDataset
import numpy as np

torch.manual_seed(42)

mean = (0.4914, 0.4822, 0.4465)
std  = (0.2023, 0.1994, 0.2010)

class SmallCNN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Conv2d(3, 32, 3, padding=1), nn.ReLU(), nn.MaxPool2d(2),
            nn.Conv2d(32, 64, 3, padding=1), nn.ReLU(), nn.MaxPool2d(2),
            nn.Flatten(),
            nn.Linear(64 * 8 * 8, 128), nn.ReLU(),
            nn.Linear(128, 10),
        )
    def forward(self, x):
        return self.net(x)

train_ds = torchvision.datasets.CIFAR10(
    root="/tmp/cifar", train=True, download=True,
    transform=T.Compose([T.ToTensor(), T.Normalize(mean, std)]))
test_ds  = torchvision.datasets.CIFAR10(
    root="/tmp/cifar", train=False, download=True,
    transform=T.Compose([T.ToTensor(), T.Normalize(mean, std)]))

train_loader = DataLoader(train_ds, batch_size=256, shuffle=True, num_workers=0)

model = SmallCNN()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
for _ in range(4):
    model.train()
    for x, y in train_loader:
        optimizer.zero_grad()
        nn.CrossEntropyLoss()(model(x), y).backward()
        optimizer.step()

def evaluate_with_noise(model, dataset, sigma):
    loader = DataLoader(dataset, batch_size=512, shuffle=False, num_workers=0)
    correct = total = 0
    model.eval()
    with torch.no_grad():
        for x, y in loader:
            x_noisy = x + sigma * torch.randn_like(x)
            preds = model(x_noisy).argmax(1)
            correct += (preds == y).sum().item()
            total += len(y)
    return 100 * correct / total

sigmas = [0.0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5]
print(f"{'Sigma':>8}  {'Accuracy (%)':>14}")
print("-" * 26)
for sigma in sigmas:
    acc = evaluate_with_noise(model, test_ds, sigma)
    bar = "#" * int(acc / 2)
    print(f"{sigma:>8.2f}  {acc:>12.1f}%  {bar}")

The output of this experiment typically shows clean accuracy near 70% for a small four-epoch model, falling to near chance (10%) by \(\sigma = 0.5\). The degradation curve is approximately sigmoidal: mild noise (\(\sigma \leq 0.1\)) has small effect, but accuracy collapses rapidly between \(\sigma = 0.15\) and \(\sigma = 0.3\). This fragility is not a consequence of model size; it reflects the high-frequency structure of natural images and the susceptibility of learned representations to perturbations that are not in the support of the training distribution.

259.7 7. Domain Generalization

Domain generalization (DG) formalizes the problem of learning from multiple source domains \(\mathcal{D}_\text{src} = \{P_1, \ldots, P_K\}\) and generalizing to an unseen target domain \(P_\text{tgt}\), with no access to target-domain data at training time. Unlike domain adaptation, which can use unlabeled target data, DG must rely entirely on the structure of the source domains.

Several approaches have been proposed. Domain randomization augments the visual domain space by applying strong random appearance transforms (color, texture, geometry), making the source distribution so broad that any realistic test distribution falls within it. This was originally developed for sim-to-real transfer in robotics (Tobin et al. 2017) but applies equally to vision tasks.

JigSaw puzzles (Carlucci et al. 2019) add a self-supervised auxiliary task of predicting the permutation of shuffled patches, encouraging the model to learn structural rather than texture features. Adaptive Data Augmentation (ADA) in DG settings learns augmentation policies that maximize domain diversity in the source distribution.

Group distributionally robust optimization (Group DRO, Sagawa et al. 2020) replaces the ERM objective with a minimax formulation:

\[\min_\theta \max_{g \in \mathcal{G}} \mathcal{L}_g(\theta)\]

where \(\mathcal{G}\) is a set of groups (subpopulations) and \(\mathcal{L}_g\) is the loss on group \(g\). Group DRO ensures that the model does not improve average loss by degrading performance on minority groups. It requires group annotations at training time, which limits its applicability when group membership is unknown.

259.8 8. Invariant Risk Minimization

Arjovsky et al. (2019) proposed Invariant Risk Minimization (IRM) as a principled framework for learning representations that are causally predictive rather than spuriously correlated. The key insight is that if a representation \(\phi: \mathcal{X} \to \mathcal{H}\) is truly invariant, then the optimal linear classifier on top of \(\phi\) should be the same across all training environments.

Formally, IRM seeks \(\phi\) and a classifier \(w\) such that \(w\) is simultaneously optimal for every environment:

\[w \in \argmin_{\bar{w}} \mathcal{L}_e(\bar{w} \circ \phi) \quad \forall e \in \mathcal{E}_\text{train}\]

This bi-level constraint is computationally intractable in its exact form. The practical relaxation penalizes the gradient of the loss with respect to a fixed scalar \(w = 1\):

\[\mathcal{R}_\text{IRM}(\phi) = \sum_{e \in \mathcal{E}_\text{train}} \left[ \mathcal{L}_e(\phi) + \lambda \left\| \nabla_{w \mid w=1} \mathcal{L}_e(w \cdot \phi) \right\|^2 \right]\]

The penalty term \(\|\nabla_{w \mid w=1} \mathcal{L}_e\|^2\) is zero when \(w=1\) is already optimal for environment \(e\), and positive otherwise. Minimizing it across environments forces the representation to make \(w=1\) approximately optimal everywhere, which is possible only when the representation has removed features whose correlations with the label differ across environments.

The theoretical guarantee (Arjovsky et al. 2019, Theorem 1) states that under linear structural causal models, IRM recovers the invariant causal features in the limit of sufficient data and environments. In practice, IRM has shown mixed empirical results: it outperforms ERM on synthetic benchmarks with known spurious correlations (Colored MNIST) but yields modest improvements on WILDS datasets. This gap between theory and practice has motivated several extensions, including IRMv1 variants, risk extrapolation (REx), and out-of-distribution generalization via causality (Peters et al. 2016).

The following diagram illustrates the relationship between spurious and invariant features under environment shift:

graph LR
    A["Input X"] --> B["Invariant features Z_I"]
    A --> C["Spurious features Z_S"]
    B --> D["Label Y"]
    E["Environment e"] --> C
    E --> F["Spurious correlation changes"]
    C --> F
    B --> G["Stable causal path"]
    G --> D

259.9 9. The Accuracy-Robustness Relationship

An empirical regularity observed across many benchmarks is that clean accuracy and corruption robustness are positively correlated across model families: larger models trained on more data tend to be more robust. Miller et al. (2021) showed that when models are evaluated on both ImageNet and ImageNet-v2 (a replicated test set with slight distribution shift), accuracy on both sets follows a tight linear relationship across dozens of architectures and training procedures.

This suggests that some of the benefit attributed to robustness-specific training methods may be confounded with overall model quality. Taori et al. (2020) found that standard augmentation and training improvements that raise clean ImageNet accuracy tend to raise ImageNet-C accuracy by a corresponding amount, placing models along a single effective robustness curve. Methods that achieve gains above this baseline curve – genuinely improving corruption robustness beyond what clean accuracy improvement alone predicts – are AugMix, DeepAugment+AugMix, and data augmentation with style transfer.

The practical implication is two-fold: first, investing in data quality, scale, and architecture is broadly beneficial for natural robustness; second, targeted robustness methods provide an additional increment of improvement over this baseline, and their evaluation must be controlled for the baseline accuracy level to avoid confounding.

259.10 10. Subpopulation Shift and Fairness

A special case of natural distribution shift is subpopulation shift: the test distribution has different group proportions than the training distribution, and the model’s accuracy varies across groups. If group \(g\) is underrepresented in training and the model relies on spurious features correlated with group membership, accuracy on \(g\) at test time (where its proportion is higher) will be systematically lower.

This is not merely a robustness failure; it is also an equity failure. Koh et al. (2021) measure both average accuracy and worst-group accuracy in WILDS, and show that ERM often achieves good average accuracy by sacrificing accuracy on minority groups. Group DRO directly addresses this by explicitly optimizing worst-group accuracy. The challenge is that group annotations are often unavailable; methods such as JTT (Liu et al. 2021) and SUBG (Zhang et al. 2022) attempt to infer group structure without explicit labels.

259.11 References

  1. Hendrycks, D. and Dietterich, T. (2019). Benchmarking Neural Network Robustness to Common Corruptions and Perturbations. International Conference on Learning Representations. arXiv:1903.12261. https://arxiv.org/abs/1903.12261

  2. Hendrycks, D., Basart, S., Mu, N., Kadavath, S., Wang, F., Dorundo, E., Desai, R., Zhu, T., Parajuli, S., Guo, M., Song, D., Steinhardt, J., and Gilmer, J. (2021). The Many Faces of Robustness: A Critical Analysis of Out-of-Distribution Generalization. International Conference on Computer Vision. arXiv:2006.16241. https://arxiv.org/abs/2006.16241

  3. Koh, P. W., Sagawa, S., Marklund, H., Xie, S. M., Zhang, M., Balsubramani, A., Hu, W., Yasunaga, M., Phillips, R. L., Gao, I., Lee, T., David, E., Stavness, I., Guo, W., Earnshaw, B., Haque, I., Beery, S. M., Leskovec, J., Kundaje, A., Pierson, E., Levine, S., Finn, C., and Liang, P. (2021). WILDS: A Benchmark of in-the-Wild Distribution Shifts. Proceedings of the 38th International Conference on Machine Learning. arXiv:2012.07421. https://arxiv.org/abs/2012.07421

  4. Hendrycks, D., Mu, N., Basart, S., Steinhardt, J., and Song, D. (2020). AugMix: A Simple Data Processing Method to Improve Robustness and Uncertainty. International Conference on Learning Representations. arXiv:1912.02781. https://arxiv.org/abs/1912.02781

  5. Arjovsky, M., Bottou, L., Gulrajani, I., and Lopez-Paz, D. (2019). Invariant Risk Minimization. arXiv:1907.02893. https://arxiv.org/abs/1907.02893

  6. Cubuk, E. D., Zoph, B., Shlens, J., and Le, Q. V. (2020). RandAugment: Practical Automated Data Augmentation with a Reduced Search Space. Advances in Neural Information Processing Systems 33. arXiv:1909.13719. https://arxiv.org/abs/1909.13719

  7. Sagawa, S., Koh, P. W., Hashimoto, T. B., and Liang, P. (2020). Distributionally Robust Neural Networks. International Conference on Learning Representations. arXiv:1911.08731. https://arxiv.org/abs/1911.08731

  8. Geirhos, R., Rubisch, P., Michaelis, C., Bethge, M., Wichmann, F. A., and Brendel, W. (2019). ImageNet-Trained CNNs Are Biased Towards Texture; Increasing Shape Bias Improves Accuracy and Robustness. International Conference on Learning Representations. arXiv:1811.12231. https://arxiv.org/abs/1811.12231

  9. Taori, R., Dave, A., Shankar, V., Carlini, N., Recht, B., and Schmidt, L. (2020). Measuring Robustness to Natural Distribution Shifts in Image Classification. Advances in Neural Information Processing Systems 33. arXiv:2007.00644. https://arxiv.org/abs/2007.00644

  10. Peters, J., Buhlmann, P., and Meinshausen, N. (2016). Causal Inference by Using Invariant Prediction: Identification and Confidence Intervals. Journal of the Royal Statistical Society: Series B, 78(5), 947–1012. https://doi.org/10.1111/rssb.12167