GLM: 稳健线性回归

GLM: 稳健线性回归

本教程是关于贝叶斯广义线性模型（GLMs）三部分系列中的第二部分，最初发表在Thomas Wiecki的博客上：

1. 线性回归

2. 稳健线性回归

3. 分层线性回归

在这篇博客文章中，我将写到：

• 少数异常值如何显著影响线性回归模型的拟合。

• 如何用学生T分布替换普通似然函数来产生稳健回归。

在线性回归教程中，我描述了最小化回归线的平方距离与最大化来自回归线的均值的正态分布的似然是相同的。这种后者的概率表达式使我们能够轻松地制定一个贝叶斯线性回归模型。

这在模拟数据上效果非常好。然而，模拟数据的问题在于，它是模拟的。在现实世界中，情况往往会变得更加复杂，像正态性这样的假设很容易被几个异常值所破坏。

让我们看看如果在上一篇文章的模拟数据中添加一些异常值会发生什么。

首先，让我们导入我们的模块。

%matplotlib inline

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor
import pytensor.tensor as pt
import xarray as xr

RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")

创建一些示例数据，但也添加一些异常值。

size = 100
true_intercept = 1
true_slope = 2

x = np.linspace(0, 1, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * x
# add noise
y = true_regression_line + rng.normal(scale=0.5, size=size)

# Add outliers
x_out = np.append(x, [0.1, 0.15, 0.2])
y_out = np.append(y, [8, 6, 9])

data = pd.DataFrame(dict(x=x_out, y=y_out))

绘制数据以及真实的回归线（左上角的三个点是我们添加的异常值）。

fig = plt.figure(figsize=(7, 5))
ax = fig.add_subplot(111, xlabel="x", ylabel="y", title="Generated data and underlying model")
ax.plot(x_out, y_out, "x", label="sampled data")
ax.plot(x, true_regression_line, label="true regression line", lw=2.0)
plt.legend(loc=0);

稳健回归

正态似然

让我们看看如果我们通过拟合一个具有正态似然的回归模型来估计我们的贝叶斯线性回归模型会发生什么。 请注意，bambi库提供了一个易于使用的工具，使得可以使用一行代码构建等效模型。 使用Bambi的此笔记本的版本可在bambi的文档中找到

with pm.Model() as model:
    xdata = pm.ConstantData("x", x_out, dims="obs_id")

    # define priors
    intercept = pm.Normal("intercept", mu=0, sigma=1)
    slope = pm.Normal("slope", mu=0, sigma=1)
    sigma = pm.HalfCauchy("sigma", beta=10)

    mu = pm.Deterministic("mu", intercept + slope * xdata, dims="obs_id")

    # define likelihood
    likelihood = pm.Normal("y", mu=mu, sigma=sigma, observed=y_out, dims="obs_id")

    # inference
    trace = pm.sample(tune=2000)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [intercept, slope, sigma]

100.00% [6000/6000 00:06<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 2_000 tune and 1_000 draw iterations (4_000 + 2_000 draws total) took 17 seconds.

为了评估拟合效果，下面的代码通过从后验分布中获取回归参数，计算后验预测回归线，并绘制其中20条回归线。

post = az.extract(trace, num_samples=20)
x_plot = xr.DataArray(np.linspace(x_out.min(), x_out.max(), 100), dims="plot_id")
lines = post["intercept"] + post["slope"] * x_plot

plt.scatter(x_out, y_out, label="data")
plt.plot(x_plot, lines.transpose(), alpha=0.4, color="C1")
plt.plot(x, true_regression_line, label="True regression line", lw=3.0, c="C2")
plt.legend(loc=0)
plt.title("Posterior predictive for normal likelihood");

如你所见，拟合结果相当偏斜，并且由于不同的后验预测回归线范围较广，我们在估计中存在相当大的不确定性。为什么会这样？原因是正态分布在尾部没有太多的质量，因此一个异常值会对拟合结果产生强烈的影响。

频率主义者会估计一个稳健回归，并使用非二次距离度量来评估拟合度。

但是，贝叶斯学派该怎么办呢？由于问题在于正态分布的尾部较轻，我们可以假设我们的数据不是正态分布的，而是根据学生T分布分布的，它具有更重的尾部，如下所示[Kruschke, 2014, Gelman 等, 2013]。

让我们看看这两个分布，以便对它们有一个直观的了解。

normal_dist = pm.Normal.dist(mu=0, sigma=1)
t_dist = pm.StudentT.dist(mu=0, lam=1, nu=1)
x_eval = np.linspace(-8, 8, 300)
plt.plot(x_eval, pm.math.exp(pm.logp(normal_dist, x_eval)).eval(), label="Normal", lw=2.0)
plt.plot(x_eval, pm.math.exp(pm.logp(t_dist, x_eval)).eval(), label="Student T", lw=2.0)
plt.xlabel("x")
plt.ylabel("Probability density")
plt.legend();

如你所见，在均值（本例中为0）远离的值的概率在T分布下比在正态分布下更可能出现。

下面是一个PyMC模型，其中likelihood项遵循StudentT分布，自由度为\(\nu=3\)，而不是Normal分布。

with pm.Model() as robust_model:
    xdata = pm.ConstantData("x", x_out, dims="obs_id")

    # define priors
    intercept = pm.Normal("intercept", mu=0, sigma=1)
    slope = pm.Normal("slope", mu=0, sigma=1)
    sigma = pm.HalfCauchy("sigma", beta=10)

    mu = pm.Deterministic("mu", intercept + slope * xdata, dims="obs_id")

    # define likelihood
    likelihood = pm.StudentT("y", mu=mu, sigma=sigma, nu=3, observed=y_out, dims="obs_id")

    # inference
    robust_trace = pm.sample(tune=4000)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [intercept, slope, sigma]

100.00% [10000/10000 00:16<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 4_000 tune and 1_000 draw iterations (8_000 + 2_000 draws total) took 24 seconds.

robust_post = az.extract(robust_trace, num_samples=20)
x_plot = xr.DataArray(np.linspace(x_out.min(), x_out.max(), 100), dims="plot_id")
robust_lines = robust_post["intercept"] + robust_post["slope"] * x_plot

plt.scatter(x_out, y_out, label="data")
plt.plot(x_plot, robust_lines.transpose(), alpha=0.4, color="C1")
plt.plot(x, true_regression_line, label="True regression line", lw=3.0, c="C2")
plt.legend(loc=0)
plt.title("Posterior predictive for Student-T likelihood")

Text(0.5, 1.0, 'Posterior predictive for Student-T likelihood')

在那里，好多了！异常值几乎不会影响我们的估计，因为我们的似然函数假设异常值比在正态分布下更有可能出现。

概述

• 通过将似然函数从正态分布改为学生T分布——后者在尾部有更多的质量——我们可以执行稳健回归。

扩展：

• 学生T分布除了均值和方差之外，还有一个称为自由度的第三个参数，用于描述尾部应放置多少质量。这里它被设置为1，这使得尾部具有最大的质量（将其设置为无穷大将得到一个正态分布！）。可以很容易地在这个参数上放置一个先验，而不是固定它，我将其留给读者作为练习 ;).

• T 分布也可以用作先验分布。参见 多层次建模的贝叶斯方法入门

• 我们如何测试我们的数据是否是正态分布的，或者在重要方面是否违反了这一假设？查看这篇精彩的博客文章，Probably Overthinking It，作者是Allen Downey。

作者

• 改编自Thomas Wiecki的博客

• 由 @fonnesbeck 于 2016 年 9 月更新（pymc#1378）

• 由 @chiral-carbon 于 2021 年 8 月更新（pymc-examples#205）

• 由Conor Hassan、Igor Kuvychko、Reshama Shaikh和Oriol Abril Pla于2022年更新

• 使用 PyMC v5 重新运行，作者：Reshama Shaikh，2023年1月

参考资料

[1]

约翰·克鲁施克。贝叶斯数据分析：R、JAGS和Stan的教程。学术出版社，2014年。

[2]

安德鲁·格尔曼, 约翰·B·卡林, 哈尔·S·斯特恩, 大卫·B·邓森, 阿基·维塔里, 和唐纳德·B·鲁宾. 贝叶斯数据分析. 查普曼和霍尔/CRC, 2013.

水印

%load_ext watermark
%watermark -n -u -v -iv -w -p xarray

Last updated: Tue Jan 10 2023

Python implementation: CPython
Python version       : 3.11.0
IPython version      : 8.8.0

xarray: 2022.12.0

matplotlib: 3.6.2
pandas    : 1.5.2
numpy     : 1.24.1
pymc      : 5.0.1
pytensor  : 2.8.11
arviz     : 0.14.0
xarray    : 2022.12.0

Watermark: 2.3.1