贝叶斯分析之利用线性回归模型理解并预测数据（三）

这一节主要讲多元线性回归模型

一元线性回归讨论的是一个因变量与一个自变量的关系，但是在很多例子中，模型可能包含多个自变量。在一元线性回归模型中，我们希望一条直线来解释数据，而在多元线性回归模型中，我们希望找到一个维度为 m 的超平面。

以 y=α+0.9∗x1+1.5∗x2y = \alpha + 0.9 * x_1 + 1.5 * x_2y=α+0.9∗x1+1.5∗x2 为例讲解多元线性回归模型
先生成数据（导入的库与第一节基本一样）

np.random.seed(314)
N = 100
alpha_real = 2.5
beta_real = [0.9, 1.5]
eps_real = np.random.normal(0, 0.5, size=N)X = np.array([np.random.normal(i, j, N) for i,j in zip([10, 2], [1, 1.5])])
X_mean = X.mean(axis=1, keepdims=True)
X_centered = X - X_mean
y = alpha_real + np.dot(beta_real, X) + eps_real

def scatter_plot(x, y):plt.figure(figsize=(10, 10))for idx, x_i in enumerate(x):plt.subplot(2, 2, idx+1)plt.scatter(x_i, y)plt.xlabel('$x_{}$'.format(idx), fontsize=16)plt.ylabel('$y$', rotation=0, fontsize=16)plt.subplot(2, 2, idx+2)plt.scatter(x[0], x[1])plt.xlabel('$x_{}$'.format(idx-1), fontsize=16)plt.ylabel('$x_{}$'.format(idx), rotation=0, fontsize=16)scatter_plot(X_centered, y)
plt.savefig('B04958_04_25.png', dpi=300, figsize=(5.5, 5.5))

再建立多元线性回归模型

with pm.Model() as model_mlr:alpha_tmp = pm.Normal('alpha_tmp', mu=0, sd=10)beta = pm.Normal('beta', mu=0, sd=1, shape=2)epsilon = pm.HalfCauchy('epsilon', 5)mu = alpha_tmp + pm.math.dot(beta, X_centered)alpha = pm.Deterministic('alpha', alpha_tmp - pm.math.dot(beta, X_mean)) y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y)start = pm.find_MAP()step = pm.NUTS(scaling=start)trace_mlr = pm.sample(5000, step=step, start=start, nchains=1)varnames = ['alpha', 'beta','epsilon']
pm.traceplot(trace_mlr, varnames)
plt.savefig('B04958_04_26.png', dpi=300, figsize=(5.5, 5.5));

这里使用 pm.math.dot() 来定义变量 mu，即线性代数中的点乘（或矩阵相乘）。

现在看一下推断出来的参数的总结

学了怎么多的线性回归模型，模型显然已经不是重点，接下来讲解多元回归模型中的三个重点：

混淆变量
多重共线性或相关性太高
隐藏的有效变量

1. 混淆变量

一个变量 z 与预测变量 x,y 都相关，如果去掉 z 之后能用 x 预测出 y，则称 z 为混淆变量。

看上面的概念阐述可能很难理解，我们从数据上来理解一下

np.random.seed(314)
N = 100
x_1 = np.random.normal(size=N)
# 其实是 x_1 决定着 x2，并且 x_1 决定着 y，所以 x_1 在分析过程中容易被忽略（混淆变量）
x_2 = x_1 + np.random.normal(size=N, scale=1)
y = x_1 + np.random.normal(size=N)
X = np.vstack((x_1, x_2))scatter_plot(X, y)
plt.savefig('B04958_04_27.png', dpi=300, figsize=(5.5, 5.5));

建立多元线性回归模型，对各个系数进行求解

可以看到 β1\beta_1β1 接近 0，这意味着 x2x_2x2 对 yyy 来说几乎没有作用（即多余的变量）。

2. 多重共线性或相关性太高

修改上述代码给 x1x_1x1 增加一个很小的扰动，因而两个变量可以看做是一样的，即 x1x_1x1 与 x2x_2x2 之间的关系是一条斜率接近于 1 的直线。

x_2 = x_1 + np.random.normal(size=N, scale=0.01)

根据 β\betaβ 系数画出2D的核密度估计图

3. 隐藏的有效变量

每个单独变量 x 不足以预测 y，如果将 x 组合在一起后就可以预测 y。因变量之间具有相关性，每个因变量都有反作用，因而忽略其中任何一个都会造成对变量影响力的低估>

先生成模拟数据，注意观察 x_0 与 x_1 的联系

np.random.seed(314)
N = 100
r = 0.8
x_0 = np.random.normal(size=N)
x_1 = np.random.normal(loc=x_0 * r, scale=(1 - r ** 2) ** 0.5)
y = np.random.normal(loc=x_0 - x_1)
X = np.vstack((x_0, x_1))scatter_plot(X, y)
plt.savefig('B04958_04_31.png', dpi=300, figsize=(5.5, 5.5));

再建立模型并求解

with pm.Model() as model_ma:alpha = pm.Normal('alpha', mu=0, sd=10)beta = pm.Normal('beta', mu=0, sd=10, shape=2)#beta = pm.Normal('beta', mu=0, sd=10)epsilon = pm.HalfCauchy('epsilon', 5)mu = alpha + pm.math.dot(beta, X)#mu = alpha + beta * X[0]y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y)start = pm.find_MAP()step = pm.NUTS(scaling=start)trace_ma = pm.sample(5000, step=step, start=start, nchains=1)pm.traceplot(trace_ma)
plt.savefig('B04958_04_32.png', dpi=300, figsize=(5.5, 5.5));

我们根据结果发现，beta 的值接近 1 和 -1，即 x1x_1x1 与 yyy 正相关，而 x2x_2x2 与 yyy 负相关。

项目源码：https://github.com/dhuQChen/BayesianAnalysis