PyMC3 贝叶斯线性回归预测与 sklearn.datasets
Posted
技术标签:
【中文标题】PyMC3 贝叶斯线性回归预测与 sklearn.datasets【英文标题】:PyMC3 Bayesian Linear Regression prediction with sklearn.datasets 【发布时间】:2016-09-15 17:16:15 【问题描述】:我一直在尝试使用 PyMC3
和来自 @ 中的数据集的 REAL DATA(即不是来自线性函数 + 高斯噪声)来实现 贝叶斯线性回归模型987654326@。我选择了形状为(442, 10)
的属性数量最少的回归数据集(即load_diabetes()
);即442 samples
和10 attributes
。
我相信我的模型工作正常,后验结果看起来足够好,可以尝试预测以弄清楚这些东西是如何工作的,但是......我意识到我不知道如何使用这些贝叶斯模型进行预测!我试图避免使用glm
和patsy
表示法,因为我很难理解使用它时实际发生了什么。
我尝试了以下操作: Generating predictions from inferred parameters in pymc3 还有http://pymc-devs.github.io/pymc3/posterior_predictive/,但我的模型要么在预测方面非常糟糕,要么我做错了。
如果我确实正确地进行了预测(我可能不是),那么任何人都可以帮助我优化我的模型。我不知道至少mean squared error
、absolute error
或类似的东西是否适用于贝叶斯框架。理想情况下,我想得到一个 number_of_rows 数组 = 我的 X_te
属性/数据测试集中的行数,以及作为后验分布样本的列数。
import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from scipy import stats, optimize
from sklearn.datasets import load_diabetes
from sklearn.cross_validation import train_test_split
from theano import shared
np.random.seed(9)
%matplotlib inline
#Load the Data
diabetes_data = load_diabetes()
X, y_ = diabetes_data.data, diabetes_data.target
#Split Data
X_tr, X_te, y_tr, y_te = train_test_split(X,y_,test_size=0.25, random_state=0)
#Shapes
X.shape, y_.shape, X_tr.shape, X_te.shape
#((442, 10), (442,), (331, 10), (111, 10))
#Preprocess data for Modeling
shA_X = shared(X_tr)
#Generate Model
linear_model = pm.Model()
with linear_model:
# Priors for unknown model parameters
alpha = pm.Normal("alpha", mu=0,sd=10)
betas = pm.Normal("betas", mu=0,#X_tr.mean(),
sd=10,
shape=X.shape[1])
sigma = pm.HalfNormal("sigma", sd=1)
# Expected value of outcome
mu = alpha + np.array([betas[j]*shA_X[:,j] for j in range(X.shape[1])]).sum()
# Likelihood (sampling distribution of observations)
likelihood = pm.Normal("likelihood", mu=mu, sd=sigma, observed=y_tr)
# Obtain starting values via Maximum A Posteriori Estimate
map_estimate = pm.find_MAP(model=linear_model, fmin=optimize.fmin_powell)
# Instantiate Sampler
step = pm.NUTS(scaling=map_estimate)
# MCMC
trace = pm.sample(1000, step, start=map_estimate, progressbar=True, njobs=1)
#Traceplot
pm.traceplot(trace)
# Prediction
shA_X.set_value(X_te)
ppc = pm.sample_ppc(trace, model=linear_model, samples=1000)
#What's the shape of this?
list(ppc.items())[0][1].shape #(1000, 111) it looks like 1000 posterior samples for the 111 test samples (X_te) I gave it
#Looks like I need to transpose it to get `X_te` samples on rows and posterior distribution samples on cols
for idx in [0,1,2,3,4,5]:
predicted_yi = list(ppc.items())[0][1].T[idx].mean()
actual_yi = y_te[idx]
print(predicted_yi, actual_yi)
# 158.646772735 321.0
# 160.054730647 215.0
# 149.457889418 127.0
# 139.875149489 64.0
# 146.75090354 175.0
# 156.124314452 275.0
【问题讨论】:
听起来不错,我完全理解。我现在就脱掉它 已经完成了,谢谢! 【参考方案1】:我认为您的模型的问题之一是您的数据具有非常不同的比例,您的“Xs”范围约为 0.3,“Ys”范围约为 300。因此,您应该期望您的先验指定的斜率(和 sigma)更大。一个合乎逻辑的选择是调整您的先验,如下例所示。
#Generate Model
linear_model = pm.Model()
with linear_model:
# Priors for unknown model parameters
alpha = pm.Normal("alpha", mu=y_tr.mean(),sd=10)
betas = pm.Normal("betas", mu=0, sd=1000, shape=X.shape[1])
sigma = pm.HalfNormal("sigma", sd=100) # you could also try with a HalfCauchy that has longer/fatter tails
mu = alpha + pm.dot(betas, X_tr.T)
likelihood = pm.Normal("likelihood", mu=mu, sd=sigma, observed=y_tr)
step = pm.NUTS()
trace = pm.sample(1000, step)
chain = trace[100:]
pm.traceplot(chain);
后验预测检查表明您有一个或多或少合理的模型。
sns.kdeplot(y_tr, alpha=0.5, lw=4, c='b')
for i in range(100):
sns.kdeplot(ppc['likelihood'][i], alpha=0.1, c='g')
另一种选择是通过标准化将数据置于相同的比例,这样做你会得到斜率应该在 +-1 左右,通常你可以对任何数据使用相同的扩散先验(有用的东西,除非你有可以使用的信息先验)。事实上,很多人都推荐这种做法用于广义线性模型。您可以在本书doing bayesian data analysis 或Statistical Rethinking 中阅读更多相关信息
如果您想预测值,您有多种选择,一种是使用推断参数的平均值,例如:
alpha_pred = chain['alpha'].mean()
betas_pred = chain['betas'].mean(axis=0)
y_pred = alpha_pred + np.dot(betas_pred, X_tr.T)
另一种选择是使用pm.sample_ppc
来获取预测值的样本,这些预测值会考虑到您估计中的不确定性。
执行 PPC 的主要思想是将预测值与您的数据进行比较,以检查它们在哪里一致,哪里不一致。例如,此信息可用于改进模型。正在做
pm.sample_ppc(trace, model=linear_model, samples=100)
将为您提供 100 个样本,每个样本具有 331 个预测观察值(因为在您的示例中,y_tr
的长度为 331)。因此,您可以将每个预测数据点与取自后验的大小为 100 的样本进行比较。您会得到预测值的分布,因为后验本身就是可能参数的分布(分布反映了不确定性)。
关于sample_ppc
的参数:samples
指定从后验得到多少个点,每个点都是参数向量。
size
指定使用该参数向量对预测值进行采样的次数(默认为 size=1
)。
您在这个tutorial 中有更多使用sample_ppc
的示例
【讨论】:
我对如何解释 sample_ppc 输出有点困惑。pm.sample_ppc(trace, model=linear_model, samples=1000)
每个 dict 元素的形状是 (1000, 111)
是我给它的 111 个测试样本(X_te)的 1000 个后验样本吗?即每个样本有 1000 个可能的预测?
samples
和size
有什么区别?【参考方案2】:
标准化 (X - u) / σ,您的自变量也可能工作得很好,因为您的 beta 方差对于所有变量都是一致的,但它们的尺度不同。
另一点可能是,如果您使用pm.math.dot
,考虑到 f(x) = intercept + Xβ + ε,计算矩阵向量积可能会更有效。
【讨论】:
以上是关于PyMC3 贝叶斯线性回归预测与 sklearn.datasets的主要内容,如果未能解决你的问题,请参考以下文章
使用 PyMC3 和大型数据集进行贝叶斯线性回归 - 括号嵌套级别超过最大值且性能缓慢