首页 文章

使用sklearn.datasets进行PyMC3贝叶斯线性回归预测

提问于
浏览 次
13

我一直在尝试使用来自 sklearn.datasets 中数据集的真实数据 PyMC3 来实现 Bayesian Linear Regression 模型（即不是来自线性函数高斯噪声） . 我选择了具有最小数量属性的回归数据集（即 load_diabetes() ），其形状为 (442, 10) ;也就是 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)

enter image description here

# 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
python statistics probability bayesian pymc3

2 回答

  • 0

    我认为你的模型存在的问题之一就是你的数据有着非常不同的比例，你的“Xs”有~0.3范围，你的“Ys”有~300 . 因此，您应该期望您的先验者指定更大的斜率（和西格玛） . 一个合理的选择是调整您的先验，如下例所示 .

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

    enter image description here

    后验预测检查表明您有一个或多或少合理的模型 .

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

    enter image description here

    另一个选择是通过标准化将数据放在相同的比例中，这样做你会得到斜率应该在-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 的示例

    回复于
  • 8

    标准化（X - u）/σ，你的自变量也可以很好地工作，因为你的beta的方差对于所有变量都是一致的，但它们有不同的比例 .

    另一点可能是，如果你使用 pm.math.dot ，考虑到f（x）=截距Xβε，计算矩阵向量积可能更有效 .

    回复于

相关问题