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

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

编辑 | sunlei

发布 | ATYUN订阅号

如果你认为贝叶斯定理是反直觉的,那么建立在贝叶斯定理基础上的贝叶斯统计就很难理解。在这一点上我和你的感受完全一致。

学习贝叶斯统计有无数的理由,尤其是贝叶斯统计正在成为表达和理解下一代深度神经网络的强大框架。

我相信,对于我们必须学习的东西,在我们能使用它们之前,我们通过使用它们来学习。生活中没有什么是如此艰难,通过我们采取的方式我们可以让它变得更容易。

所以,这是我简化它的方法:与其在开始时使用过多的理论或术语,不如让我们关注贝叶斯分析的机制,特别是如何使用PyMC3和ArviZ进行贝叶斯分析和可视化。在记忆无穷无尽的术语之前,我们将对解决方案进行编码并将结果可视化,并使用术语和理论解释模型。

PyMC3是一个用于概率编程的Python库,语法非常简单直观。ArviZ是一个与PyMC3携手工作的Python库,它可以帮助我们解释和可视化后验分布。

我们将把贝叶斯方法应用到一个实际问题中,展示一个端到端的贝叶斯分析,它从构建问题到建立模型到获得先验概率再到在Python中实现最终的后验分布。

在我们开始之前,让我们先得出一些基本的直觉:

贝叶斯模型也被称为概率模型,因为它们是用概率建立的。贝叶斯利用概率作为量化不确定性的工具。因此,我们得到的答案是分布而不是点估计。

贝叶斯方法步骤

步骤1:建立关于数据的信念,包括先验函数和似然函数。

步骤2:根据我们对数据的信念,使用数据和概率,更新我们的模型,检查我们的模型是否与原始数据一致。

步骤3:根据模型更新数据视图。

数据

由于我对使用机器学习进行价格优化很感兴趣,所以我决定将贝叶斯方法应用到一个西班牙高铁票价数据集,该数据集可以在这里找到。感谢Gurus团队对数据集的搜集。

from scipy import stats
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import seaborn as sns
import pandas as pd
from theano import shared
from sklearn import preprocessingprint('Running on PyMC3 v{}'.format(pm.__version__))data = pd.read_csv('renfe.csv')
data.drop('Unnamed: 0', axis = 1, inplace=True)
data = data.sample(frac=0.01, random_state=99)
data.head(3)
data.isnull().sum()/len(data)

价格栏中有12%的值丢失了,我决定用相应票价类型的平均值来填充它们。还用最常见的值填充其他两个分类列。

data['train_class'] = data['train_class'].fillna(data['train_class'].mode().iloc[0])
data['fare'] = data['fare'].fillna(data['fare'].mode().iloc[0])
data['price'] = data.groupby('fare').transform(lambda x:    

x.fillna(x.mean()))

高斯推论

az.plot_kde(data['price'].values, rug=True)
plt.yticks([0], alpha=0);

铁路票价的KDE图显示了高斯分布,除了几十个远离均值的数据点。

假设高斯分布是对火车票价格的恰当描述。由于我们不知道平均值或标准差,所以必须为这两个值设置优先级。因此,一个合理的模型可以是这样的。

模型

我们将对票价数据进行高斯推断。这里有一些模型选择。

我们将在PyMC3中这样实例化模型:

  • PyMC3中的模型规范封装在with语句中。

先验选择:

  • μ,指人口。正态分布很广。我不知道μ的可能的值,我可以设置先验。根据经验,我知道火车票价格不能低于0或高于300,所以我将均匀分布的边界设为0和300。你可能有不同的经历,设定不同的界限。完全没问题。如果你有比我更可靠的信息,请使用它!
  • σ,标准差的人口。只能是正的,因此使用半正态分布。再来一次,非常宽广。

票价似然函数的选择:

  • y是一个观测变量,代表的数据来自正态分布的参数μ、σ。
  • 使用螺母取样绘制1000个后验样本。

使用PyMC3,我们可以将模型写成:

with pm.Model() as model_g:
    μ = pm.Uniform('μ', lower=0, upper=300)
    σ = pm.HalfNormal('σ', sd=10)
    y = pm.Normal('y', mu=μ, sd=σ, observed=data['price'].values)
    trace_g = pm.sample(1000, tune=1000)

y表示可能性。这就是我们告诉PyMC3我们要根据已知(数据)为未知条件设置条件的方式。

我们绘制高斯模型轨迹。这是运行在一个Theano图表下的引擎盖。

az.plot_trace(trace_g);
  • 在左边,我们有一个KDE图,对于x轴上的每个参数值我们在y轴上得到一个概率它告诉我们参数值的可能性有多大。
  • 在右边,我们得到了采样过程中每个步骤的单独采样值。从轨迹图中,我们可以从后面直观地得到可信的值。
  • 上面的图中每个参数都有一行。对于这个模型,后面是二维的,因此上图显示了每个参数的边缘分布。

这里有几点需要注意:

  • 我们对单个参数的采样链(左)似乎很好地收敛和稳定(没有大的漂移或其他奇怪的模式)。
  • 每个变量的最大后验估计(左侧分布的峰值)非常接近真实参数。

我们可以绘制参数的联合分布。

az.plot_joint(trace_g, kind='kde', fill_last=False);

我看不出这两个参数之间有任何关联。这意味着模型中可能没有共线性。这是很好的。

我们还可以对每个参数的后验分布进行详细的总结。

az.summary(trace_g)

我们还可以通过生成一个分布的均值和最高后验密度(HPD)的图来直观地看到上述总结,并解释和报告贝叶斯推断的结果。

az.plot_posterior(trace_g);
  • 与频域推理不同,在贝叶斯推理中,我们得到了整个值的分布。
  • 每次ArviZ计算和报告HPD时,默认情况下它将使用94%的值。
  • 请注意,HPD间隔与confidence间隔不同。
  • 在这里,我们可以这样解释,94%的概率,相信平均票价在8欧元到64.4欧元之间。

我们可以用Gelman Rubin检验来验证链的收敛性。接近1.0的值表示收敛。

pm.gelman_rubin(trace_g)
bfmi = pm.bfmi(trace_g)
max_gr = max(np.max(gr_stats) for gr_stats in pm.gelman_rubin(trace_g).values())
(pm.energyplot(trace_g, legend=False, figsize=(6, 4)).set_title("BFMI = {}\nGelman-Rubin = {}".format(bfmi, max_gr)));

我们的模型收敛得很好,Gelman Rubin数据看起来也不错。

在今天的学习当中,我们了解了贝叶斯方法步骤和高斯推论,也将贝叶斯方法应用到一个实际问题中,展示一个端到端的贝叶斯分析,明天我会继续更新接下来的内容。

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

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

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

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 数据科学家常犯的十大编程错误

    数据科学家是“比任何软件工程师都更擅长统计,比任何软件工程师都更擅长软件工程的的统计学家”。许多数据科学家都有统计学背景却缺乏在软件工程方面的经验。我是资深的数...

    AiTechYun
  • 如何规避线性回归的陷阱(下)

    在上一部分中,我们学习了线性回归的概念和规避线性回归陷阱的前两个解决方案,今天我们继续学习剩余的两个方案。

    AiTechYun
  • 【学术】在C ++中使用TensorFlow训练深度神经网络

    你可能知道TensorFlow的核心是用C++构建的,然而只有python的API才能获得多种便利。 当我写上一篇文章时,目标是仅使用TensorFlow的C ...

    AiTechYun
  • 手把手 | 教你爬下100部电影数据:R语言网页爬取入门指南

    大数据文摘
  • Python中的数据标准化

    数据标准化 数据标准化是指将数据按比例缩放,使之落入到特定区间。 为了消除量纲的影响,方便进行不同变量间的比较分析。 0-1标准化: x=(x-min)/(ma...

    Erin
  • jquery mobile 移动web(2)

    button 按钮   data-role="button" 将超链接变成button。   具有icon 图标的button 组件...

    用户1197315
  • 寻找和为定值的两个数

    输入一个整数数组和一个整数,在数组中查找一对数,满足他们的和正好是输入的那个整数,如果有多对数的和等于输入的整数,则全部输出,要求输出的结果中不应该出现重复,如...

    陌无崖
  • 12种用于Python数据分析的Pandas技巧

    本文将介绍12种用于数据分析的Pandas技巧,为了更好地描述它们的效果,这里我们用一个数据集辅助进行操作。

    崔庆才
  • 深入浅出谈「大数据」| MTdata小讲堂

    美图数据技术团队
  • R分类算法-神经网络算法

    神经网络(Artifical Neural Network) 神经网络(人工神经网络),是一种模仿生物网络(动物的中枢神经系统,特别是大脑)的结构和功能的数学模...

    Erin

扫码关注云+社区

领取腾讯云代金券