分层部分池化

假设你被要求估计几位棒球选手的击球技能。其中一个这样的表现指标是击球率。由于选手们参加的比赛数量不同，并且在击球顺序中的位置也不同，每位选手的击球次数也不同。然而，你希望估计所有选手的技能，包括那些击球机会相对较少的选手。

所以，假设一个球员只上场击球4次，并且从未击中球。他们是一个糟糕的球员吗？

作为免责声明，本笔记本的作者对棒球及其规则的了解很少，甚至可以说几乎没有。他一生中击球的次数大约是“4”次。

数据

我们将使用棒球数据 [Efron 和 Morris, 1975]。

方法

我们将使用PyMC来估计每位球员的击球率。在估计了数据集中所有球员的平均击球率后，我们可以利用这些信息来估计一位新球员的击球率，该球员的数据较少（即 4次击球）。

在没有贝叶斯分层模型的情况下，这个问题有两种方法：

1. 独立计算每位球员的击球率（无池化）

2. 在假设每个人的潜在平均值相同的情况下计算总体平均值（完全池化）

当然，这两种方法都不现实。显然，所有球员的击球技能并不相同，因此全球平均值是不合理的。同时，职业棒球运动员在许多方面是相似的，因此他们的平均值也不是完全独立的。

可能可以对“相似”的玩家群体进行聚类，并估计群体平均值，但使用分层建模方法是一种自然的信息共享方式，不需要识别特别的聚类。

分层部分池化的思想是建模全局表现，并使用该估计来参数化一个球员群体，该群体考虑了球员表现之间的差异。这种全局和个体表现之间的权衡将由模型自动调整。此外，由于每个球员的击球次数（即信息）不同而导致的不确定性将通过将这些估计值向全局平均值收缩来自动考虑。

关于更深入的讨论，请参阅Stan的教程 [Carpenter 等, 2016]。该模型和参数值取自该示例。

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt

%matplotlib inline

RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")

现在我们可以使用 pandas 加载数据集：

data = pd.read_csv(pm.get_data("efron-morris-75-data.tsv"), sep="\t")
at_bats, hits = data[["At-Bats", "Hits"]].to_numpy().T

现在让我们为这些数据开发一个生成模型。

我们将假设存在一个与所有球员（不仅限于我们的18名球员）预期表现相关的隐藏因子（phi）。由于总体均值是一个介于0和1之间的未知值，因此它必须在上下限之间进行界定。此外，我们假设对全球平均水平一无所知。因此，一个自然的选择是先验分布为均匀分布。

接下来，我们引入一个超参数 kappa 来考虑总体击球率的方差，我们将使用一个有界的帕累托分布。这将确保估计值落在合理的范围内。这些超参数将依次用于参数化一个beta分布，该分布非常适合对单位区间上的量进行建模。beta分布通常通过尺度参数和形状参数进行参数化，它也可以通过其均值 \(\mu \in [0,1]\) 和样本大小（方差的代理）\(\nu = \alpha + \beta (\nu > 0)\) 进行参数化。

最后一步是为每个玩家的数据（命中或未命中）指定一个采样分布，使用二项分布。这是数据应用于模型的部分。

我们可以使用 pm.Pareto('kappa', m=1.5) 来定义我们对 kappa 的先验，但帕累托分布的尾部非常长。正确探索这些尾部对于采样器来说很困难，因此我们使用指数分布的等效但更快的参数化方法。我们利用了帕累托分布的随机变量的对数遵循指数分布这一事实。

N = len(hits)
player_names = data["FirstName"] + " " + data["LastName"]
coords = {"player_names": player_names.tolist()}

with pm.Model(coords=coords) as baseball_model:
    phi = pm.Uniform("phi", lower=0.0, upper=1.0)

    kappa_log = pm.Exponential("kappa_log", lam=1.5)
    kappa = pm.Deterministic("kappa", pt.exp(kappa_log))

    theta = pm.Beta("theta", alpha=phi * kappa, beta=(1.0 - phi) * kappa, dims="player_names")
    y = pm.Binomial("y", n=at_bats, p=theta, dims="player_names", observed=hits)

回想一下，我们最初的问题是关于一个只有4次击球机会且没有安打的球员的真实击球率。我们可以将此作为模型中的一个附加变量

with baseball_model:
    theta_new = pm.Beta("theta_new", alpha=phi * kappa, beta=(1.0 - phi) * kappa)
    y_new = pm.Binomial("y_new", n=4, p=theta_new, observed=0)

该模型可以这样可视化

pm.model_to_graphviz(baseball_model)

我们现在可以使用MCMC拟合模型：

with baseball_model:
    idata = pm.sample(2000, tune=2000, chains=2, target_accept=0.95)

    # check convergence diagnostics
    assert all(az.rhat(idata) < 1.03)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [phi, kappa_log, theta, theta_new]

100.00% [8000/8000 00:26<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 2_000 tune and 2_000 draw iterations (4_000 + 4_000 draws total) took 34 seconds.

现在我们可以绘制参数的后验分布。首先，是总体超参数：

az.plot_trace(idata, var_names=["phi", "kappa"]);

因此，总体平均击球率在0.22-0.31范围内，预期值约为0.26。

接下来，数据集中所有18名球员的估计值：

az.plot_forest(idata, var_names="theta");

最后，让我们为我们的0投0中球员获取估计值：

az.plot_trace(idata, var_names=["theta_new"]);

请注意，尽管我们的额外球员没有获得任何安打，但他平均值的估计并不是零——零甚至不是一个高度可能的值。这是因为我们假设该球员是从一个由我们估计的超参数指定的分布的总体球员中抽取的。然而，该球员的估计平均值处于我们数据集中球员平均值的较低端，这表明这4次打数为估计提供了一些信息。

作者

• 由Vladislavs Dovgalecs于2016年11月撰写（pymc#1546）

• 由Adrian Seybolt于2017年6月更新（pymc#2288）

• 由Christian Luhmann于2020年8月更新（pymc#4068）

• 使用 PyMC v5 由 Reshama Shaikh 于 2023 年 1 月运行

参考资料

[1]

布拉德利·埃夫隆和卡尔·莫里斯。使用斯坦因估计量及其推广的数据分析。美国统计协会杂志，70(350):311–319, 1975。

[2]

Bob Carpenter, J Gabry, 和 B Goodrich. 分层部分池化用于重复二元试验。2016年。

水印

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

Last updated: Sat Jan 28 2023

Python implementation: CPython
Python version       : 3.11.0
IPython version      : 8.8.0

xarray: 2022.12.0

matplotlib: 3.6.2
pandas    : 1.5.2
arviz     : 0.14.0
pytensor  : 2.8.11
pymc      : 5.0.1
numpy     : 1.24.1

Watermark: 2.3.1

许可证声明

本示例库中的所有笔记本均在MIT许可证下提供，该许可证允许修改和重新分发，前提是保留版权和许可证声明。

引用 PyMC 示例

要引用此笔记本，请使用Zenodo为pymc-examples仓库提供的DOI。

重要

许多笔记本是从其他来源改编的：博客、书籍……在这种情况下，您应该引用原始来源。

同时记得引用代码中使用的相关库。

这是一个BibTeX的引用模板：

@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}

渲染后可能看起来像：