Press "Enter" to skip to content

Pymc简介及用于描述统计模型的语言

由Paul Nicklen拍摄,来自Treehugger

在我们之前关于为什么大多数贝叶斯推理的例子都对其进行了错误描述的文章中,我们澄清了贝叶斯统计学初学者之间的一个常见误解。即,贝叶斯统计学并不是通过使用贝叶斯定理来定义的,而是通过使用概率分布来表征不确定性并考虑可能结果的全范围来定义的。因此,例如,我们可以考虑这样的情况,即不是被告知在您已经确定自己已经感染的情况下,某个医疗设备在检测疾病方面的有效率为95%(我们也可以称之为真阳性率(TPR)),而是考虑该设备在检测疾病方面的有效率只有24%,69%或91%以及这些数字如何改变您在从设备上阳性测试后感染该疾病的概率。通过拒绝接受点估计(如静态TPR),而是倾向于一种概率方法,通过考虑一系列结果,我们更加忠实于贝叶斯概率的观点。

说到这里,我们在学习过程中,为了达到对概率论有贝叶斯思维方式的状态,接下来的阶段是构建整合概率方法来表征不确定性的贝叶斯统计模型。其中一个最好的工具就是专门用于进行贝叶斯推理和构建概率机器学习模型的PyMC Python包。

然而,在我们开始构建贝叶斯模型之前,我们需要进行一个快速的绕道,并介绍从业者用于构建这些模型的语言(即数学符号)。理解构建模型的语言对于两个原因是实用的:

原因一。第一个原因是学习从业者(如科学家和研究人员)如何使用应用统计模型符号来构建模型与我们在PyMC上编写可执行代码的方式直接相关。使用模型符号允许统计建模者用几行简短的符号来传达他们对模型的基本假设。我们认为,与其让观众记住复杂的词汇(如同方差性)以了解他们正在使用的模型的特性,不如使用符号更简单的沟通形式。

原因二。另一个原因是这种语言涵盖了所有领域,从艺术史到天文学和保护生物学(McElreath,2020)。因此,通过熟悉描述模型的语言,我们将在科学界中扩展我们理解和沟通知识的边界(McElreath,2020)。

作为对应用统计模型符号的一般介绍,McElreath教授(2020)将一些原则描述为科学术语描述和理解模型的基础:

  1. 从业者识别一组要处理的变量,有时这些变量是可观测的。这些被称为数据。
    a. 数据产生的不可观测现象,如平均值,被定义为参数。
  2. 从业者将每个变量定义为其他变量的函数或概率分布。
  3. 变量及其概率分布的组合定义了一个联合生成模型,我们可以用它来模拟一个假设的观察结果并分析真实的观察结果。

现在我们已经了解了模型符号世界背后的原理,让我们用一个玩具例子来说明我们如何使用它们。假设您是一位对在野外遇到温哥华岛沿海狼(Canis lupus crassodon)中的雄性狼的概率进行量化感兴趣的保育人士。

顺便一提,我们对在示例中使用温哥华岛沿海狼的原因之一是因为它们最近被发现并被归类为一个独立的亚种。因此,与其他狼亚种相比,它们的数据和信息量有限,因此在这个亚种的特征以及未来的文章中生成估计值将是有趣的。它们的数据不仅有限,而且在我们看来,它们是非常酷的动物,与其他狼亚种不同的是,它们能够游泳长距离并将海鲜纳入其天然饮食(Muñoz‐Fuentes等,2009)。从更广泛的视角来看,动物生态学和种群统计研究对动物的保护状况、种群管理和对其在某个地区更广泛的生态系统中的角色的理解都具有影响。

回到我们的例子,仅仅基于这个简单的研究问题,我们现在可以生成假设并建立一个模型,量化温哥华岛沿海狼的性别比例。首先,我们可以假设这个思想实验的结果的分布是二项分布,这意味着作为我们研究问题的结果,只能有两种可能的结果发生。要么我们在穿越温哥华岛的旅程中遇到一只雄性沿海狼,要么是一只雌性。实际上,我们可以使用以下代码模拟当两种二项结果都有50%的机会发生时,我们期望的二项分布的样子:

# 导入所需的库
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import scipy.stats as stats
az.style.use("arviz-darkgrid")

x_binom = np.random.binomial(n=10, p=0.5, size=100)
sns.histplot(x_binom, kde=True, discrete=True)
plt.xlabel("硬币正面朝上的次数(N=10)")
plt.ylabel("样本数量")
plt.title('可视化二项分布')

Pymc简介及用于描述统计模型的语言 四海 第2张

从图中,我们可以找到一个接近正态分布的近似结果,这是等同于进行了10次抛硬币实验的结果,x轴描述了硬币正面朝上的次数。如果我们像上面的代码一样重复这个实验100次,y轴详细说明了我们在x轴上得到的结果出现的频率。换句话说,在我们重复进行这个实验的100次中,大约有30次实验的结果是在10次抛硬币中大约有5次正面朝上。请注意,每次执行上面的代码时,图的结果会稍微有所变化。

让我们在我们的模型中加入另一个假设,我们假设任何一种结果发生的机会都是相等的。在二项分布中,并不总是这样,所以这不是一个我们应该默认采取的假设。在分布中,当任何和所有结果发生的机会都是相等的时,被称为均匀分布。它通常被用作模型符号,以传达对模型中可能结果的零先验知识。然而,随着足够的实践和经验的积累,我们会发现,建立模型的贝叶斯从业者很少会对他们打算建模的现象没有任何先验知识。如果您好奇均匀分布在图上是什么样子,我们也可以模拟这个分布:

x_uni = np.linspace(-.25, 1.25, 100)
plt.plot(x_uni, stats.uniform.pdf(x_uni, 0, 1))
plt.xlabel("先验概率(性别比例)")
plt.ylabel("比例的可信度")
plt.title("可视化均匀分布")

Pymc简介及用于描述统计模型的语言 四海 第3张

一旦我们确定了沿海狼性别比例模型的假设,让我们在应用统计模型符号中将其正式化:

Pymc简介及用于描述统计模型的语言 四海 第4张

  • M ~ 二项分布(N, p):我们可以将第一行解释为“我们将遇到的雄性数量(M)服从二项分布,样本大小为N,先验概率为p。”这行代码代表我们的似然函数,也被称为贝叶斯定理中的证据的概率P(E | H),它是观察到给定结果的可能性。在这种情况下,“给定结果”将是我们对沿海狼的样本数据的结果,我们将很快详细介绍。
  • p ~ 均匀分布(0, 1):最后,我们可以将第二行解释为“先验概率 P(H) 在0到1之间均匀分布。”

请注意~符号表示“分布为”,并表示左侧模型的参数与其分布具有随机关系(即“随机确定”或具有“随机概率分布”属性,可以在统计上进行预测但不是精确的)。换句话说,每当我们看到~符号表示随机关系时,我们应该期望右侧数学陈述的输出是由左侧变量(在我们的情况下,要么是p要么是M)捕获的分布。表达随机关系与我们习惯于看到的数学陈述形成了鲜明对比,单个数字(而不是分布)是由等号=符号组成的方程的输出,例如我们经常在线性回归中看到的一个常见方程,即y = mx + b

在我们将模型转化为代码之前,我们需要生成一个观察样本,用于我们的似然函数(P(E | H)),它表示我们观察到给定事件结果的可能性。例如,如果我们使用Muñoz-Fuentes和她的同事(2009年)在温哥华岛沿海狼研究中获得的样本数据,我们会发现在28只带有标记性别的狼中,有15只是雄性。在这种情况下,我们可以将我们的似然函数描述为在我们的样本数据集中遇到15只雄性狼的概率。

在确定了我们的变量之后,我们现在可以使用PyMC直接将我们刚刚形式化的模型符号化为代码:

# 从Muñoz-Fuentes等人(2009年)的研究中读取我们的样本数据
wolf_samples = "https://raw.githubusercontent.com/vanislekahuna/Statistical-Rethinking-PyMC/main/Data/Vancouver-Island-wolf-samples.csv"
islandwolves_df = pd.read_csv(wolf_samples)
islandwolves_df

# 确定数据集中的雄性与雌性比例
print(f"温哥华岛沿海狼性别的值计数为:")
print(islandwolves_df["Sexb"].value_counts()

# 生成我们模型的样本观察
males = islandwolves_df["Sexb"].value_counts()["M"]
females = islandwolves_df["Sexb"].value_counts()["F"]
total_samples = males + females

### 我们的PyMC模型 ###
with pm.Model() as wolf_ratio:
    # 先验值的可能结果的均匀分布
    priors = pm.Uniform("uniform_prior", lower=0, upper=1)
    # 从我们的先验均匀分布和样本数据中生成二项分布,其中28只狼中有15只是雄性
    n_males = pm.Binomial("sex_ratio", p=priors, observed=males, n=total_samples)
    # 从模型中抽取样本生成我们的后验分布
    trace_wolves = pm.sample(1000, tune=1000)
######################

现在我们的模型已经实现了,让我们不要陷入细节中,忘记构建模型的唯一目的。我们构建模型的原因是为了通过可能性分布来量化遇到雄性狼的概率,并根据其可能性对其进行排名。换句话说,我们对观察到的28只狼中的15只雄性的概率(我们的似然函数)感兴趣,给定每种可行的雄性与雌性比例,即先验分布(例如,野外所有狼中1%为雄性,2%,3%,以此类推直到100%)。这里的最终结果应该是一个显示基于我们观察到的样本数据,哪些雄性与雌性比例最有可能的后验概率分布。从上面的代码中,后验分布被捕获在`trace_wolves`对象中,它通过对我们构建的模型进行采样来生成结果的后验概率。如果你想知道,用于生成后验分布的采样方法是一个我们可以在以后的文章中讨论的主题,因为PyMC提供了几种选择,每种选择都是一个深入且广泛的兔子洞。

下面左图可视化了我们的后验分布的结果,其中x轴代表性别比例的先验概率分布,y轴代表该比例的可能性:

az.plot_trace(trace_wolves)

Pymc简介及用于描述统计模型的语言 四海 第5张

如果将刚刚构建的`wolf_ratio`模型结构化为贝叶斯定理中的变量,我们可以这样做:

Pymc简介及用于描述统计模型的语言 四海 第6张

再次强调,这里的目标是生成每个先验值在其分布中的可能性的后验分布。我们并不试图找到具有最高概率的性别比例。话虽如此,如果你对该值感兴趣,我们可以使用PyMC的`find_map`函数来给出最有可能的先验值(即最有可能的雄性与雌性比例)在后验分布中:

pm.find_MAP(model=wolf_ratio)

或者使用Arviz的`summary`函数,该函数不仅提供后验分布的最大先验概率(MAP),还提供包含89%分布的最高密度区间(有时也称为兼容区间),这两个先验值包含了分布的89%。

az.summary(trace_wolves, round_to=2, kind="stats", hdi_prob=0.89)

以上就是我们对大多数学术研究中使用的统计模型符号语言进行基本介绍的内容。然后我们将这些模型符号整合到我们的第一个PyMC模型中,通过该模型我们可以模拟在野外遇到一只温哥华岛沿海雄性狼的可能性,假设我们对它们的雄性和雌性性别比一无所知。熟悉模型符号将为使用PyMC构建更复杂的贝叶斯模型奠定坚实的基础。我们计划在不久的将来发布更多关于这个主题的博客文章,首先是关于贝叶斯线性预测的文章,敬请关注!

参考文献:

McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with examples in R and Stan. Routledge.

Muñoz‐Fuentes, V., Darimont, C. T., Wayne, R. K., Paquet, P. C., & Leonard, J. A. (2009). Ecological factors drive differentiation in wolves from British Columbia. Journal of Biogeography, 36(8), 1516–1531. https://www.raincoast.org/files/publications/papers/Munoz_et_al_2009_JBiogeog.pdf

Muñoz-Fuentes, V., Darimont, C.T., Paquet, P.C., & Leonard, J.A. (2010). The genetic legacy of extirpation and re-colonization in Vancouver Island wolves. Conservation Genetics, 11, 547–556. https://digital.csic.es/bitstream/10261/58619/1/evol.pdf

Leave a Reply

Your email address will not be published. Required fields are marked *