首页 > 编程语言 >Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据

时间:2023-11-01 11:34:13浏览次数:48  
标签:BMA 后验 Python 模型 贝叶斯 PyMC3 使用 我们 pm

当面对多个模型时,我们有多种选择。模型选择因其简单性而具有吸引力,但我们正在丢弃有关模型中不确定性的信息。

 

 

print(f"Runing

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_元模型

模型平均

一种替代方法是执行模型选择,但讨论所有不同的模型以及给定信息准则的计算值。重要的是要将所有这些数字和测试放在我们问题的背景下,以便我们和我们的客户能够更好地了解我们方法可能存在的局限性和缺点。如果你在学术界,你可以使用这种方法向论文、演示文稿、论文等的讨论部分添加元素。

另一种方法是执行模型平均。现在的想法是使用模型的加权平均值生成元模型(和元预测)。有几种方法可以做到这一点,PyMC3 包括其中的 3 种,我们将简要讨论,您将在 Yuling Yao 等人的工作中找到更彻底的解释。

伪贝叶斯模型平均

贝叶斯模型可以通过其边缘概率进行加权,这被称为贝叶斯模型平均。我们可以使用以下公式来做到这一点:

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_数据_02

这种方法称为伪贝叶斯模型平均或类似赤池的加权,是一种启发式方法,用于根据信息标准值计算每个模型(给定一组固定的模型)的相对概率。看看分母只是一个归一化项,以确保权重总和为 1。

使用贝叶斯自举进行伪贝叶斯模型平均

上述计算权重的公式是一种非常好且简单的方法,但有一个主要警告,它没有考虑 IC 计算中的不确定性。

堆叠

在PyMC3中实现的第三种方法被称为预测分布的堆叠,并且最近被提出。我们希望在一个元模型中组合多个模型,以最小化元模型和真实生成模型之间的分歧,当使用对数评分规则时,这相当于:

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_数据_03

加权后验预测样本

一旦我们计算了权重,使用上述 3 种方法中的任何一种,我们就可以使用它们来获得加权后验预测样本。PyMC3 提供了以简单方式执行这些步骤的函数,因此让我们通过示例查看它们的实际效果。

简而言之,我们的问题如下:我们想探索几种灵长类动物的乳汁成分,假设来自大脑较大的灵长类动物的雌性产生更有营养的牛奶(这样做是为了*支持这种大大脑的发育)。对于进化生物学家来说,这是一个重要的问题,为了给出和回答,我们将使用3个变量,两个预测变量:新皮层的比例与总质量的比较 大脑和母亲体重的对数。对于预测变量,每克牛奶的千卡。使用这些变量,我们将构建 3 个不同的线性模型:

  1. 仅使用新皮层变量的模型
  2. 仅使用质量变量对数的模型
  3. 使用两个变量的模型

 

 

d.iloc[:, 1:] = d.iloc[:, 1:] - d.iloc[:, 1:].mean()
d.head()

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_权重_04

现在我们有了数据,我们将仅使用 neocortex

 

 

with pm.Model() as model_0:
  
    trace_0 = pm.sample(2000, return_inferencedata=True)

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_元模型_05

第二个模型与第一个模型完全相同,只是我们现在使用质量的对数

 

 

with pm.Model() as model_1:
 

    trace_1 = pm.sample(2000, return_inferencedata=True)

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_数据_06

最后是第三个模型使用 neocortex和 变量log_mass

 

 

with pm.Model() as model_2:
   

    trace_2 = pm.sample(2000, return_inferencedata=True)

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_权重_07

现在我们已经对 3 个模型的后验进行了采样,我们将对它们进行视觉比较。一种选择是使用forestplot支持绘制多个迹线的函数。

 

 

az.plot_fo

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_数据_08

 另一种选择是在同一图中绘制多条迹线是使用densityplot 。

 

 

az.plot_d

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_元模型_09

 现在我们已经对 3 个模型的后验进行了采样,我们将使用 WAIC(广泛适用的信息标准)来比较 3 个模型。我们可以使用 PyMC3 附带的compare功能来做到这一点。

 

 

comp = az.compare(model_dict)
comp

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_数据_10

我们可以看到最好的模型是,具有两个预测变量的模型。请注意,数据帧按从最低到最高 WAIC 的顺序(最差的模型)。

现在,我们将使用copmuted来生成预测,而不是基于单个模型,而是基于加权模型集。

 

 

ppc_w = pm.sample_posterior_predictive_w(

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_数据_11

请注意,我们正在传递按其索引排序的权重。

我们还将计算最低 WAIC 模型的 PPC

 

 

ppc_2 = pm.sample_posterior_predi

比较这两种预测的一种简单方法是绘制它们的平均值和 hpd 区间

 

 

plt.yticks([])
plt.ylim(-1, 2)
plt.legend();

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_元模型_12

正如我们所看到的,两个预测的平均值几乎相同,但加权模型中的不确定性更大。我们已经有效地将我们应该选择哪个模型的不确定性传递到后验预测样本中。

结语:

还有其他方法可以平均模型,例如,显式构建一个包含我们拥有的所有模型的元模型。然后,我们在模型之间跳转时执行参数推理。这种方法的一个问题是,在模型之间跳跃可能会妨碍后验的正确采样。

版本信息

 

 

%load_ext watermark
%watermark -n -u -v -iv -w

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_元模型_13

Python用PyMC3贝叶斯模型平均BMA:采样、信息准则比较和预测可视化灵长类动物的乳汁成分数据_权重_14



标签:BMA,后验,Python,模型,贝叶斯,PyMC3,使用,我们,pm
From: https://blog.51cto.com/u_14293657/8120293

相关文章

  • Python用PyMC贝叶斯GLM广义线性模型、NUTS采样器拟合、后验分布可视化
    尽管贝叶斯方法相对于频率主义方法的理论优势已经在其他地方进行了详细讨论,但其更广泛采用的主要障碍是“可用性”。而使用贝叶斯方法,客户可以按照自己认为合适的方式定义模型。线性回归在此示例中,我们将帮助客户从最简单的GLM–线性回归开始。一般来说,频率论者对线性回归的看......
  • 21.13 Python 实现端口流量转发
    端口流量转发(PortForwarding)是一种网络通信技术,用于将特定的网络流量从一个端口或网络地址转发到另一个端口或地址。它在网络中扮演着一个非常重要的角色,在Python语言中实现端口转发非常容易。如下这段代码实现了一个基本的TCP端口映射,将本地指定端口的流量转发到指定的远程IP和......
  • Python学习笔记(一)蒙特卡罗算法求圆周率π
    绪论\(\pi\)(圆周率)是数学和物理学普遍存在的常数之一,可以被定义为圆周长和直径之比或者圆的面积与半径平方之比(\(l=2\pir\)和\(S=\pir^2\))。\(\pi\)是一个无理数,下面将用蒙特卡罗算法求\(\pi\)的数值近似。要求1.要求能算到小数点后面越多越好‪‬‪‬‪‬‪‬‪‬‮‬‪‬‫......
  • 【python基础】repr函数
     描述repr()函数将对象转化为供解释器读取的形式。语法以下是repr()方法的语法:repr(object)参数object--对象。返回值返回一个对象的string格式。实例#coding=UTF-8s="物品\t单价\t数量\n包子\t1\t2"print(s)print(repr(s))output:物品单价数量包子1......
  • 阿基米德螺线和花瓣曲线的Python绘制
    frommatplotlibimportpyplotaspltimportnumpyasnpdefplt_draw_dynamic(x,y):foriinrange(x.shape[0]):plt.pause(0.001)plt.cla()plt.xlim((-200,200))plt.ylim((-200,200))plt.plot(x[:i],y[:i])......
  • python深浅拷贝学习
    copy的原文链接(仅供自己学习查看):python浅析格式化输出和深浅copy-战争热诚-博客园(cnblogs.com) 首先我们从切片技术说起。它应用于所有的序列,包括:列表,字符串,元祖。但是切片不能应用于字典,对于字典只能使用D.copy()和D.deepcopy()方法。下面具体说一下深......
  • 带有最小间隔时间的队列读取实现 —— Python异步编程
     (注:照片源自免费网站,地址:https://www.freepik.com/photos/angry-panda/13)  ==================================================   ==================================================......
  • Python求Π
    fromrandomimportrandomfrommathimportsqrtfromtimeimportperf_counterimporttime#进度条print("2022310143117")foriinrange(1,101):print("\r",end="")print("进度:{}%:".format(i),"▓"*(i//......
  • 如何在linux系统中安装python3.8.1 并卸载 python3.6.2 更新python3引导到3.8.1
    安装python3.8.1步骤1:检查Python版本在终端中输入以下命令来检查当前安装的Python版本:python--version步骤2:安装编译Python所需的依赖项更新系统软件包,并安装构建Python所需的一些工具和库。在终端中运行以下命令:sudoaptupdatesudoapt-getinstall-ybuild-essen......
  • 【爬虫实战】用Python采集任意小红书笔记下的评论,爬了10000多条,含二级评论!
    目录一、爬取目标二、爬虫代码讲解2.1分析过程2.2爬虫代码三、演示视频一、爬取目标您好!我是@马哥python说,一名10年程序猿。我们继续分享Python爬虫的案例,今天爬取小红书上指定笔记("巴勒斯坦"相关笔记)下的评论数据。老规矩,先展示结果:截图1:截图2:截图3:共爬取了1w多条"......