python 贝叶斯回归模型 贝叶斯回归预测

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

问题：利用给定的数据建立 y 与 x 之间的线性模型 $y=\alpha+\beta x$

1. 构造出数据集

先导入相应的一系列库

%matplotlib inline
import pymc3 as pm
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
palette = 'muted'
sns.set_palette(palette); sns.set_color_codes(palette)
np.set_printoptions(precision=2)
pd.set_option('display.precision', 2)
sns.set()

假设一个线性模型： $y=2.5+ 0.9 * x$，在生成数据时加一个扰动项（eps_real）

np.random.seed(1)
N = 100
alfa_real = 2.5
beta_real = 0.9
eps_real = np.random.normal(0, 0.5, size=N)

x = np.random.normal(10, 1, N)
y_real = alfa_real + beta_real * x 
y = y_real + eps_real

numpy.random.normal(loc=0.0, scale=1.0, size=None)
loc：float, 此概率分布的均值（对应着整个分布的中心center）
scale：float，此概率分布的标准差（对应于分布的宽度，scale越大越矮胖，scale越小越瘦高）
size：int or tuple of ints，输出的shape，默认为None，只输出一个值（注意是个长度为 size 的向量）

现在看看数据分布与真实的线性回归模型

plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(x, y, 'b.')
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.plot(x, y_real, 'k')
plt.subplot(1,2,2)
sns.kdeplot(y)
plt.xlabel('$y$', fontsize=16)
plt.tight_layout()
plt.savefig('B04958_04_02.png', dpi=300, figsize=(5.5, 5.5))

2. 建立线性回归模型

with pm.Model() as model:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=1)
    epsilon = pm.HalfCauchy('epsilon', 5)

    mu = pm.Deterministic('mu', alpha + beta * x)
    y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y)
    
    start = pm.find_MAP() 
    step = pm.Metropolis() 
    trace = pm.sample(10000, step, start, nchains=1)
pm.traceplot(trace)

先看看模型拟合的效果，左图为核密度估计（Kernel Density Estimation，KDE）图，右图描述了每一步采样过程中得到的采样值。

注意：

1. $\epsilon$
2. KDE图做了归一化（使得概率密度积分为1），因此纵轴数值有很大差别。

再分析一下 $\alpha$, $\beta$ 自相关性：

可以发现 $\alpha$, $\beta$ 有很糟糕的自相关性。

请思考：为什么 $\alpha$, $\beta$ 因为在使用“最小二乘法”的条件下，不论用哪条拟合的直线都会经过一点，采样点的均值点$（ \bar{x},\bar{y}）$

---

基于均方误差（SE，即对误差平方求和）最小化来进行模型求解的方法称为“最小二乘法”。在线性回归中，最小二乘法就是试图找到一条直线，使得所有样本到直线上的欧式距离之和最小。
那么这个求解方法到底是什么呢？其实就是对 $w$, $b$ 求导！
$\frac{\partial E_{(w,b)}} {\partial w} =2(w\sum_{i=1}^{m}x_i^2-\sum_{i=1}^{m}(y_i-b)x_i)$$\frac{\partial E_{(w,b)}} {\partial b} =2(mb-\sum_{i=1}^{m}(y_i - wx_i))$
进行求解：
$w = \frac{ \sum_{i=1}^{m} y_i(x_i-\bar x)}{\sum_{i=1}^{m} x_i^2- \frac{1}{m}(\sum_{i=1}^{m} x_i)^2}$$b= \frac{1}{m}\sum_{i=1}^{m} (y_i - wx_i)$
将求得的 $w$, $b$ 带入模型有：$\bar y = w\bar x + b$

---

由上可知，拟合直线的过程相当于将直线固定在数据的中心（均值点）上进行旋转，斜率越大截距越小，因此两个参数是相关的。
现将后验画出来，可以发现 $\alpha$, $\beta$

sns.kdeplot(trace['alpha'], trace['beta'])
plt.xlabel(r'$\alpha$', fontsize=16)
plt.ylabel(r'$\beta$', fontsize=16, rotation=0)
plt.savefig('B04958_04_05.png', dpi=300, figsize=(5.5, 5.5));

3. 解决高自相关性

3.1 中心化或者标准化

$x^{$
中心化使得 $x^{$ 的中心在 0 附近，从而使得修改斜率时旋转点变成了截距点，参数空间也会变得不那么自相关。直观上解释就是不再必经均值点，即 $\bar y = w\bar x + b$

$x^{$$y^{$
标准化好处之一是我们对数据使用了相同的弱先验，而不必关心数据的具体值域有多大。

3.2 更换采样方法

1. NUTS算法
2. Metropolis算法

with pm.Model() as model_n:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=1)
    epsilon = pm.HalfCauchy('epsilon', 5)
    
    mu = pm.Deterministic('mu', alpha + beta * x)

    y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y)
    
    start = pm.find_MAP() 
    # 更改采样方法？
    step = pm.NUTS() 
    trace_n = pm.sample(2000, step=step, start=start, nchains=1)
pm.traceplot(trace_n, varnames)

看看模型拟合的效果，可以发现各个参数都有较好的混合度。

现在再看看参数的自相关性与最后结果

4. 对后验进行解释和可视化

注意：以下后验均是基于更改采样方法（采用NUTS算法）得出的后验。

4.1 后验的不确定性

plt.plot(x, y, 'b.');

idx = range(0, len(trace_n['alpha']), 10)
plt.plot(x, trace_n['alpha'][idx] + trace_n['beta'][idx] *  x[:,np.newaxis], c='gray', alpha=0.5);

plt.plot(x, alpha_m + beta_m * x, c='k', label='y = {:.2f} + {:.2f} * x'.format(alpha_m, beta_m))

plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.legend(loc=2, fontsize=14)
plt.show()

其中半透明的的直线表示后验的不确定性，可以发现中间部分的不确定性较低，不过直线并没有相交于一个点（后验并不强制所有的直线都穿过均值点）。

4.2 预测值$\hat y$的 HPD 区间

ppc = pm.sample_ppc(trace_n, samples=1000, model=model_n)

plt.plot(x, y, 'b.')
plt.plot(x, alpha_m + beta_m * x, c='k', label='y = {:.2f} + {:.2f} * x'.format(alpha_m, beta_m))

sig0 = pm.hpd(ppc['y_pred'], alpha=0.5)[idx]
sig1 = pm.hpd(ppc['y_pred'], alpha=0.05)[idx]
plt.fill_between(x_ord, sig0[:,0], sig0[:,1], color='gray', alpha=1)
plt.fill_between(x_ord, sig1[:,0], sig1[:,1], color='gray', alpha=0.5)

plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.show()

对预测值进行采用，将50%HPD区间用深灰色区域表示，将95%HPD区间用浅灰色表示。

5. 皮尔逊相关系数

对于皮尔逊相关系数 r，我们需要了解以下几点：

1. 衡量两个变量之间线性相关性（因此 r = 0 表示没有线性关系，可能存在其他非线性关系）；
2. $r= \beta\frac{\sigma(x)}{\sigma(y)}$，其中 $\beta$
3. 决定系数是皮尔逊相关系数的平方。

5.1 利用PyMC3计算 r

with pm.Model() as model_n:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=1)
    epsilon = pm.HalfCauchy('epsilon', 5)
    
    mu = alpha + beta * x
    y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y)

    rb = pm.Deterministic('rb', (beta * x.std() / y.std()) ** 2)

    y_mean = y.mean()
    ss_reg = pm.math.sum((mu - y_mean) ** 2)
    ss_tot = pm.math.sum((y - y_mean) ** 2)
    rss = pm.Deterministic('rss', ss_reg/ss_tot)

    start = pm.find_MAP()
    step = pm.NUTS()
    trace_n = pm.sample(2000, step=step, start=start, nchains=1)
pm.traceplot(trace_n)

5.2 根据多元高斯分布计算 r

先看一下高斯分布，其中标准差 $\sigma=1$。

sigma_x1 = 1
sigmas_x2 = [1, 2]
rhos = [-0.99, -0.5, 0, 0.5, 0.99]

k, l = np.mgrid[-5:5:.1, -5:5:.1]
pos = np.empty(k.shape + (2,))
pos[:, :, 0] = k; pos[:, :, 1] = l

f, ax = plt.subplots(len(sigmas_x2), len(rhos), sharex=True, sharey=True, figsize=(12, 8))
# f.figure(figsize=(5, 1))
for i in range(2):
    for j in range(5):
        sigma_x2 = sigmas_x2[i]
        rho = rhos[j]
        cov = [[sigma_x1**2, sigma_x1*sigma_x2*rho], [sigma_x1*sigma_x2*rho, sigma_x2**2]]
        rv = stats.multivariate_normal([0, 0], cov)
        ax[i,j].contour(k, l, rv.pdf(pos))
        ax[i,j].plot(0, 0, 
        label="$\\sigma_{{x2}}$ = {:3.2f}\n$\\rho$ = {:3.2f}".format(sigma_x2, rho),  alpha=0)
        ax[i,j].legend()
ax[1,2].set_xlabel('$x_1$')
ax[1,0].set_ylabel('$x_2$')
plt.show()

由于并不知道这些标准差（组成的协方差矩阵），因此我们可以通过设置标准差的先验，然后利用这些值手动构造协方差矩阵。

data = np.stack((x, y)).T

with pm.Model() as pearson_model:
    
    mu = pm.Normal('mu', mu=data.mean(0), sd=10, shape=2)
    
    sigma_1 = pm.HalfNormal('simga_1', 10)
    sigma_2 = pm.HalfNormal('sigma_2', 10)
    rho = pm.Uniform('rho', -1, 1)
    
    cov = pm.math.stack(([sigma_1**2, sigma_1*sigma_2*rho], [sigma_1*sigma_2*rho, sigma_2**2]))
    
    y_pred = pm.MvNormal('y_pred', mu=mu, cov=cov, observed=data)
    
    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    trace_p = pm.sample(1000, step=step, start=start, nchains=1)

pm.traceplot(trace_p)
plt.show()