贝叶斯回归R语言代码 贝叶斯回归案例分析

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

---

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

---

以 $y = \alpha + 0.9 * x_1 + 1.5 * x_2$
先生成数据（导入的库与第一节基本一样）

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. 混淆变量
2. 多重共线性或相关性太高
3. 隐藏的有效变量

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));

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

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

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

修改上述代码给 $x_1$ 增加一个很小的扰动，因而两个变量可以看做是一样的，即 $x_1$ 与 $x_2$

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，即 $x_1$ 与 $y$ 正相关，而 $x_2$ 与 $y$

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