这一节主要讲多元线性回归模型
一元线性回归讨论的是一个因变量与一个自变量的关系，但是在很多例子中，模型可能包含多个自变量。在一元线性回归模型中，我们希望一条直线来解释数据，而在多元线性回归模型中，我们希望找到一个维度为 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 = 100alpha_real = 2.5beta_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_meany = 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 = 100x_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 = 100r = 0.8x_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 负相关。
项目源码：/dhuQChen/BayesianAnalysis
