专栏首页ATYUN订阅号Python/PyMC3/ArviZ贝叶斯统计实战(下)

Python/PyMC3/ArviZ贝叶斯统计实战(下)

编辑 | sunlei

发布 | ATYUN订阅号

在上半部分中,我们了解了贝叶斯方法步骤和高斯推论,也将贝叶斯方法应用到一个实际问题中,今天我主要介绍贝叶斯在Python中实现最终的后验分布。

前文回顾:Python/PyMC3/ArviZ贝叶斯统计实战(上)

后验预测检验(PPCs)是验证模型的一种很好的方法。其思想是使用来自后验图的参数从模型中生成数据。

现在我们已经计算了后验,我们将说明如何使用模拟结果来推导预测。

下面的函数将从跟踪中随机抽取1000个参数样本。然后,对于每个样本,它将从该样本中μ和σ值指定的正态分布中提取25798个随机数。

后预测检查

function DisplayWindowSize(){

  var w=window.innerWidth
  || document.documentElement.clientWidth
  || document.body.clientWidth;

现在,ppc包含1000个生成的数据集(每个数据集包含25798个样本),每个数据集使用与后验不同的参数设置。

_, ax = plt.subplots(figsize=(10, 5))
ax.hist([y.mean() for y in ppc['y']], bins=19, alpha=0.5)
ax.axvline(data.price.mean())
ax.set(title='Posterior predictive of the mean', xlabel='mean(x)', ylabel='Frequency');

推断的平均值与实际的火车票价格平均值非常接近。

分组比较

我们可能对不同票价类型下的价格比较感兴趣。我们将着重于估计效应的大小,即量化两类票价之间的差异。为了比较票价类别,我们将使用每种票价类型的平均值。因为我们是贝叶斯,所以我们将努力获得票价类别之间的均值差异的后验分布。

我们创建了三个变量:

价格变量,表示票价。

idx变量,一个用数字编码票价类别的分类虚拟变量。

最后是组变量,包含票价类别的数量(6)

price = data['price'].values
idx = pd.Categorical(data['fare'],
                                categories=['Flexible', 'Promo', 'Promo +',

'Adulto ida', 'Mesa', 'Individual-Flexible']).codes
groups = len(np.unique(idx))

对于6个组(票价类别),很难绘制每个组的μ和σ的跟踪图。因此,我们创建一个汇总表:

flat_fares = az.from_pymc3(trace=trace_groups)
fares_gaussian = az.summary(flat_fares)
fares_gaussian

很明显,不同组别(即票价类别)的平均票价有显著差异。

为了更清楚地说明这一点,我们在不重复比较的情况下绘制了每个票价类别之间的差异。

<li><p style="text-align: left;">Cohen’s d是比较两种方法的合适的效应大小。Cohen ‘s d通过标准差来引入每一组的变异性。</p><

><li><p style="text-align: left;">优势概率(ps)定义为随机从一个组取的数据点大于随机从另一个组取的数据点的概率。</p><>

dist = stats.norm()

_, ax = plt.subplots(5, 2, figsize=(20, 12), constrained_layout=True 

comparisons = [(i, j) for i in range(6) for j in range(i+1, 6)]
pos = [(k, l) for k in range(5) for l in (0, 1)]

for (i, j), (k, l) in zip(comparisons, pos):
    means_diff = trace_groups['μ'][:, i] - trace_groups['μ'][:, j]

    d_cohen = (means_diff / np.sqrt((trace_groups['σ'][:, i]**2 + trace_groups['σ'][:, j]**2) / 2)).mean()
    ps = dist.cdf(d_cohen/(2**0.5))
    az.plot_posterior(means_diff, ref_val=0, ax=ax[k, l])
    ax[k, l].set_title(f'$\mu_{i}-\mu_{j}$')
    ax[k, l].plot(
        0, label=f"Cohen's d = {d_cohen:.2f}\nProb sup = {ps:.2f}", alpha=0)
    ax[k, l].legend();

基本上,上面的图告诉我们,在上面的比较案例中,94%的HPD都没有包含0的参考值。这意味着对于所有的例子,我们可以排除0的差。6.1欧元到63.5欧元的平均差价足够大,顾客可以根据不同的票价类别购买机票。

贝叶斯层次线性回归

我们要建立一个模型来估计每种火车类型的火车票价格,同时估计所有火车类型的价格。这种类型的模型称为层次模型或多级模型。

01
def replace_fare(fare):

    if fare == 'Adulto ida':
        return 1
    elif fare == 'Promo +':
        return 2
    elif fare == 'Promo':
        return 3
    elif fare == 'Flexible':
        return 4
    elif fare == 'Individual-Flexible':
        return 5
    elif fare == 'Mesa':
        return 6

data['fare_encode'] = data['fare'].apply(lambda x: replace_fare(x))

label_encoder = preprocessing.LabelEncoder()
data['train_type_encode']= label_encoder.fit_transform(data['train_type'])

train_type_names = data.train_type.unique()
train_type_idx = data.train_type_encode.values

n_train_types = len(data.train_type.unique())

data[['train_type', 'price', 'fare_encode']].head()

我们将建模的数据的相关部分如下所示。我们感兴趣的是不同的火车类型是否会影响票价。

层次模型

with pm.Model() as hierarchical_model:
    # global model parameters
    α_μ_tmp = pm.Normal('α_μ_tmp', mu=0., sd=100)
    α_σ_tmp = pm.HalfNormal('α_σ_tmp', 5.)
    β_μ = pm.Normal('β_μ', mu=0., sd=100)
    β_σ = pm.HalfNormal('β_σ', 5.)

    # train type specific model parameters
    α_tmp = pm.Normal('α_tmp', mu=α_μ_tmp, sd=α_σ_tmp, shape=n_train_types) 
    # Intercept for each train type, distributed around train type mean
    β = pm.Normal('β', mu=β_μ, sd=β_σ, shape=n_train_types)
    # Model error
    eps = pm.HalfCauchy('eps', 5.)

    fare_est = α_tmp[train_type_idx] + β[train_type_idx]*data.fare_encode.values

    # Data likelihood
    fare_like = pm.Normal('fare_like', mu=fare_est, sd=eps, observed=data.price)

with hierarchical_model:
    hierarchical_trace = pm.sample(2000, tune=2000, target_accept=.9)

pm.traceplot(hierarchical_trace, var_names=['α_μ_tmp', 'β_μ', 'α_σ_tmp', 'β_σ', 'eps']);

左列是高度的边缘后验信息,告诉我们“α_μ_tmp”组平均价格水平,“β_μ”告诉我们, 与票价类型“Adulto ida”和购买票价类别“Promo +”相比,购买票价类别“Promo +”显著提高了价格。与票价类型““Promo +”等(零下无质量)相比,价格明显降低。

pm.traceplot(hierarchical_trace, var_names=['α_tmp'], coords={'α_tmp_dim_0': range(5)});

在16种火车类型中,我们可能想看看5种火车类型在票价方面的比较。从“α-tmp”的边缘可以看出,列车类型之间的价格差异很大;不同的宽度与我们对每个参数估计的置信度有关-每种列车类型的测量值越多,我们的置信度就越高。

对我们的一些估计进行不确定性量化是贝叶斯模型的一个强大功能。对于不同类型的火车,我们有一个贝叶斯可信区间。

az.plot_forest(hierarchical_trace, var_names=['α_tmp', 'β'], combined=True);

最后,我们要计算r的平方:

ppc = pm.sample_posterior_predictive(hierarchical_trace,

samples=2000, model=hierarchical_model)
az.r2_score(data.price.values, ppc['fare_like'])

这篇文章的目的是学习、实践和解释贝叶斯,而不是从数据集中产生尽可能的最佳结果。否则,我们将直接使用XGBoost。

Jupyter notebook可以在Github上找到,享受这一周剩下的时间吧。

Github链接地址:

https://github.com/susanli2016/Machine-Learning-with-Python/blob/master/Bayesian%20Statistics%20Python_PyMC3_ArviZ.ipynb

本文分享自微信公众号 - ATYUN订阅号(atyun_com),作者:关注人工智能的

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-09-08

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 【业界】iOS Bug解决办法:如何防止Siri读出隐藏的通知

    AiTechYun 编辑:yuxiangyu ? 尽管苹果尚未针对iOS错误展开修复,使得Siri能够读出隐藏的锁屏通知,但以下方法可以解决目前的安全漏洞。 ...

    AiTechYun
  • 使用keras创建一个简单的生成式对抗网络(GAN)

    然而,有些恶意的顾客为了获得金钱而出售假酒。在这种情况下,店主必须能够区分假酒和正品葡萄酒。

    AiTechYun
  • 使用Apache MXNet分类交通标志图像

    有许多深度学习的框架,例如TensorFlow、Keras、Torch和Caffe,Apache MXNet由于其在多个GPU上的可伸缩性而受到欢迎。在这篇博文...

    AiTechYun
  • 腾讯WeTest的小程序兼容测试实践之路

    ? 作 者 朱永俊,腾讯IEG高级工程师 商业转载请联系腾讯WeTest获得授权,非商业转载请注明出处。 作者导读 为了提升对微信小程序的测试能力, 腾讯W...

    WeTest质量开放平台团队
  • Python爬虫之JS逆向入门篇

    我们都知道现在是大数据时代,用爬虫来获取数据的越来越多,与之对应的就是破解反爬的难度也越来越大

    Python编程与实战
  • 面试官常说,培训机构出身的程序员“代码不干净”,是什么意思?

    现在很多企业对于培训出来的程序员都带着有色眼镜在看,甚至一些过激的企业直接把培训出来的程序员排除在外,这种做法很明显是不正确的,主要很多培训机构的宣传以及包装对...

    程序员互动联盟
  • Junit单元测试遇见的一个枚举类型的坑(枚举类型详解)

    枚举类型很早就在计算机语言中存在了,主要被用来将一组相似的值包含进一种类型中,这种类型的名称被定义成独一无二的类型描述符,这就是枚举类型。

    Criss@陈磊
  • mstsc远程报:这可能是由于CredSSP 加密Oracle修正的两种完美解决方法

    查看win10系统升级日志,果然找到了原因,是因为CVE-2018-0886 的 CredSSP 2018 年 5 月 8 日更新默认设置从“易受攻击”更改为“...

    依乐祝
  • 自定义对话框绑定控件

    forrestlin
  • 使用PHP生成ICO图标

    今天教大家如何使用PHP生成ico图标,ico图标在每个网站中都需要用到的,使用方法也是很简单的,基本上以下面的方式为主,还有其他的方式。

    小白程序猿

扫码关注云+社区

领取腾讯云代金券