这一节主要讲一元线性回归模型
问题:利用给定的数据建立 y 与 x 之间的线性模型
y
=
α
+
β
x
y=\alpha+\beta x
y=α+β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 y=2.5+ 0.9 * x 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)图,右图描述了每一步采样过程中得到的采样值。
注意:
- ϵ \epsilon ϵ 具有较好的混合度;
- KDE图做了归一化(使得概率密度积分为1),因此纵轴数值有很大差别。
再分析一下
α
\alpha
α,
β
\beta
β 自相关性:
可以发现
α
\alpha
α,
β
\beta
β 有很糟糕的自相关性。
请思考:为什么
α
\alpha
α,
β
\beta
β 自相关性很强? 因为在使用“最小二乘法”的条件下,不论用哪条拟合的直线都会经过一点,采样点的均值点
(
x
ˉ
,
y
ˉ
)
( \bar{x},\bar{y})
(xˉ,yˉ)
基于均方误差(SE,即对误差平方求和)最小化来进行模型求解的方法称为“最小二乘法”。在线性回归中,最小二乘法就是试图找到一条直线,使得所有样本到直线上的欧式距离之和最小。
那么这个求解方法到底是什么呢?其实就是对
w
w
w,
b
b
b 求导!
∂
E
(
w
,
b
)
∂
w
=
2
(
w
∑
i
=
1
m
x
i
2
−
∑
i
=
1
m
(
y
i
−
b
)
x
i
)
\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)
∂w∂E(w,b)=2(wi=1∑mxi2−i=1∑m(yi−b)xi)
∂
E
(
w
,
b
)
∂
b
=
2
(
m
b
−
∑
i
=
1
m
(
y
i
−
w
x
i
)
)
\frac{\partial E_{(w,b)}} {\partial b} =2(mb-\sum_{i=1}^{m}(y_i - wx_i))
∂b∂E(w,b)=2(mb−i=1∑m(yi−wxi))
进行求解:
w
=
∑
i
=
1
m
y
i
(
x
i
−
x
ˉ
)
∑
i
=
1
m
x
i
2
−
1
m
(
∑
i
=
1
m
x
i
)
2
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}
w=∑i=1mxi2−m1(∑i=1mxi)2∑i=1myi(xi−xˉ)
b
=
1
m
∑
i
=
1
m
(
y
i
−
w
x
i
)
b= \frac{1}{m}\sum_{i=1}^{m} (y_i - wx_i)
b=m1i=1∑m(yi−wxi)
将求得的
w
w
w,
b
b
b 带入模型有:
y
ˉ
=
w
x
ˉ
+
b
\bar y = w\bar x + b
yˉ=wxˉ+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
−
x
ˉ
x^{'}=x-\bar{x}
x′=x−xˉ
中心化使得
x
′
x^{'}
x′ 的中心在 0 附近,从而使得修改斜率时旋转点变成了截距点,参数空间也会变得不那么自相关。直观上解释就是不再必经均值点,即
y
ˉ
=
w
x
ˉ
+
b
\bar y = w\bar x + b
yˉ=wxˉ+b 不再成立。
x
′
=
x
−
x
′
x
s
d
x^{'}=\frac{x-x^{'}}{x_{sd}}
x′=xsdx−x′
y
′
=
y
−
y
′
y
s
d
y^{'}=\frac{y-y^{'}}{y_{sd}}
y′=ysdy−y′
标准化好处之一是我们对数据使用了相同的弱先验,而不必关心数据的具体值域有多大。
3.2 更换采样方法
- NUTS算法
- 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 预测值 y ^ \hat y 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,我们需要了解以下几点:
- 衡量两个变量之间线性相关性(因此 r = 0 表示没有线性关系,可能存在其他非线性关系);
- r = β σ ( x ) σ ( y ) r= \beta\frac{\sigma(x)}{\sigma(y)} r=βσ(y)σ(x),其中 β \beta β 为斜率;
- 决定系数是皮尔逊相关系数的平方。
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
先看一下高斯分布,其中标准差 σ = 1 \sigma=1 σ=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()
项目源码:https://github.com/dhuQChen/BayesianAnalysis