PyMC3贝叶斯线性回归预测与sklearn.datasets

Question

我一直在尝试使用PyMC3和sklearn.datasets数据集中的真实数据（即不是线性函数+高斯噪声）来实现Bayesian Linear Regression模型。我选择了属性最少的回归数据集（即load_diabetes()），其形状为(442,10)，即442个样本和10个属性。我相信我的模型已经能够工作了，后验分布看起来足够好，可以尝试进行预测，以弄清楚这些东西的工作原理，但是......我意识到我不知道如何使用这些贝叶斯模型进行预测！我试图避免使用glm和patsy符号，因为我很难理解在使用它时实际上正在发生什么。我尝试过参考生成pymc3推断参数的预测和http://pymc-devs.github.io/pymc3/posterior_predictive/，但我的模型要么极其糟糕地进行预测，要么我做错了。如果我实际上正在正确进行预测（我很可能没有），那么有谁能帮助我优化我的模型？我不知道在贝叶斯框架中是否使用最小均方误差、绝对误差或类似的东西。理想情况下，我想获得一个数组，行数=我的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 

- O.rka

---

好的，我完全明白。我现在就去处理。 - O.rka

已经完成了，谢谢！ - halfer

Answer 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，并且通常情况下可以对任何数据使用相同的漫散先验（除非您有信息先验可用）。事实上,许多人推荐在广义线性模型中采用此方法。 您可以在书籍《贝叶斯数据分析》或《统计思维》中了解更多信息。

如果您想要预测值，则有几个选项之一是使用推断参数的平均值，例如：

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）。

在此教程中，您可以找到更多使用sample_ppc的示例。

- aloctavodia

---

我对如何解释sample_ppc输出有些困惑。pm.sample_ppc(trace, model=linear_model, samples=1000) 对于每个字典元素，形状为(1000, 111)，它是针对我提供的111个测试样本（X_te）的1000个后验样本吗？即每个样本有1000个可能的预测值？ - O.rka

“samples”和“size”的区别是什么？ - O.rka

Answer 2

标准化 (X - u) / σ，您的自变量也可能效果良好，因为所有变量的beta方差是均匀的，但它们具有不同的比例尺。

另一个要点是，如果使用 pm.math.dot，在计算矩阵向量乘积时可能更有效，因为 f(x) = intercept + Xβ + ε。

- minggli