高尔夫推杆的模型构建与扩展

这使用了并紧密遵循Andrew Gelman的案例研究，该研究是用Stan编写的。我们增加了一些新的可视化内容，并避免了使用不恰当的先验分布，但非常感谢他和Stan团队提供的精彩案例研究和软件。

我们使用来自“统计学：贝叶斯视角”的数据集 [Berry, 1996]。该数据集描述了职业高尔夫球手从不同距离推杆的结果，并且数据量足够小，我们可以直接打印并加载它，而不需要进行任何特殊的 csv 读取。

注意

本笔记本使用了不是 PyMC 依赖项的库，因此需要专门安装这些库才能运行此笔记本。打开下面的下拉菜单以获取更多指导。

Extra dependencies install instructions

为了在本地或binder上运行此笔记本，您不仅需要一个安装了所有可选依赖项的PyMC工作环境，还需要安装一些额外的依赖项。有关安装PyMC本身的建议，请参阅安装

您可以使用您喜欢的包管理器安装这些依赖项，我们提供了以下作为示例的pip和conda命令。

$ pip install

请注意，如果您想（或需要）从笔记本内部而不是命令行安装软件包，您可以通过运行pip命令的变体来安装软件包：

导入系统

!{sys.executable} -m pip install

您不应运行!pip install，因为它可能会在不同的环境中安装包，即使安装了也无法从Jupyter笔记本中使用。

另一个替代方案是使用conda：

$ conda 安装

在安装科学计算Python包时，我们推荐使用conda forge

import io

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
import scipy
import scipy.stats as st
import xarray as xr

from xarray_einstats.stats import XrContinuousRV

RANDOM_SEED = 8927
az.style.use("arviz-darkgrid")

# golf putting data from berry (1996)
golf_data = """distance tries successes
2 1443 1346
3 694 577
4 455 337
5 353 208
6 272 149
7 256 136
8 240 111
9 217 69
10 200 67
11 237 75
12 202 52
13 192 46
14 174 54
15 167 28
16 201 27
17 195 31
18 191 33
19 147 20
20 152 24"""


golf_data = pd.read_csv(io.StringIO(golf_data), sep=" ", dtype={"distance": "float"})

BALL_RADIUS = (1.68 / 2) / 12
CUP_RADIUS = (4.25 / 2) / 12

我们开始绘制数据以更好地了解其外观。隐藏单元格包含绘图代码

Show code cell source Hide code cell source

def plot_golf_data(golf_data, ax=None, color="C0"):
    """Utility function to standardize a pretty plotting of the golf data."""
    if ax is None:
        _, ax = plt.subplots()
    bg_color = ax.get_facecolor()
    rv = st.beta(golf_data.successes, golf_data.tries - golf_data.successes)
    ax.vlines(golf_data.distance, *rv.interval(0.68), label=None, color=color)
    ax.plot(
        golf_data.distance,
        golf_data.successes / golf_data.tries,
        "o",
        mec=color,
        mfc=bg_color,
        label=None,
    )

    ax.set_xlabel("Distance from hole")
    ax.set_ylabel("Percent of putts made")
    ax.set_ylim(bottom=0, top=1)

    ax.set_xlim(left=0)
    ax.grid(True, axis="y", alpha=0.7)
    return ax

ax = plot_golf_data(golf_data)
ax.set_title("Overview of data from Berry (1996)");

绘图后，我们发现一般来说，高尔夫球手在更远的地方准确性较低。请注意，这些数据是预先聚合的：我们可能能够使用逐杆数据进行更有趣的工作。此数据集似乎已被分箱到最近的英尺。

我们可能会考虑用这些数据进行预测：拟合一条曲线到这些数据上，可以让我们对中间距离做出合理的猜测，甚至可能外推到更长的距离。

Logit 模型

首先我们将拟合一个传统的logit-二项式模型。我们直接对成功的次数进行建模，使用

\[\begin{split} a, b \sim \mathcal{N}(0, 1) \\ p(\text{success}) = \operatorname{logit}^{-1}(a \cdot \text{distance} + b) \\ \text{num. successes} \sim \operatorname{Binomial}(\text{tries}, p(\text{success})) \end{split}\]

以下是如何在PyMC中编写该模型。我们在模型变量中使用下划线后缀，以避免污染命名空间。我们还使用pymc.MutableData，以便稍后当我们处理新数据集时可以替换数据。

with pm.Model() as logit_model:
    distance_ = pm.MutableData("distance", golf_data["distance"], dims="obs_id")
    tries_ = pm.MutableData("tries", golf_data["tries"], dims="obs_id")
    successes_ = pm.MutableData("successes", golf_data["successes"], dims="obs_id")

    a_ = pm.Normal("a")
    b_ = pm.Normal("b")

    pm.Binomial(
        "success",
        n=tries_,
        p=pm.math.invlogit(a_ * distance_ + b_),
        observed=successes_,
        dims="obs_id",
    )


pm.model_to_graphviz(logit_model)

我们有一些直觉认为\(a\)应该是负的，并且\(b\)应该是正的（因为当\(\text{distance} = 0\)时，我们期望能接近100%的推杆成功率）。不过，我们并没有将这些放入模型中。我们将其作为基线，并且我们也可以等待，看看是否需要添加更强的先验。

with logit_model:
    logit_trace = pm.sample(1000, tune=1000, target_accept=0.9)

az.summary(logit_trace)

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

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

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
a	-0.255	0.007	-0.268	-0.243	0.000	0.000	880.0	1256.0	1.0
b	2.223	0.057	2.114	2.324	0.002	0.001	931.0	1259.0	1.0

我们看到 \(a\) 和 \(b\) 的符号与我们预期的一致。采样器没有发出任何不良警告。查看摘要，有效样本数量是合理的，rhat 接近 1。这是一个小模型，所以我们不会过于仔细地检查拟合情况。

我们绘制了50个\(p(\text{success})\)的后验样本以及期望值。此外，我们还从后验预测中抽取了500个点进行绘制：

# Draw posterior predictive samples
with logit_model:
    logit_trace.extend(pm.sample_posterior_predictive(logit_trace))

# hard to plot more than 400 sensibly
logit_post = az.extract(logit_trace, num_samples=400)
logit_ppc = az.extract(logit_trace, group="posterior_predictive", num_samples=400)
const_data = logit_trace["constant_data"]

logit_ppc_success = logit_ppc["success"] / const_data["tries"]

# Plotting
ax = plot_golf_data(golf_data)
t_ary = np.linspace(CUP_RADIUS - BALL_RADIUS, golf_data.distance.max(), 200)
t = xr.DataArray(t_ary, coords=[("distance", t_ary)])
logit_post["expit"] = scipy.special.expit(logit_post["a"] * t + logit_post["b"])

ax.plot(
    t,
    logit_post["expit"].T,
    lw=1,
    color="C1",
    alpha=0.5,
)

ax.plot(t, logit_post["expit"].mean(dim="sample"), color="C2")

ax.plot(golf_data.distance, logit_ppc_success, "k.", alpha=0.01)
ax.set_title("Logit mean and posterior predictive");

100.00% [4000/4000 00:00<00:00]

拟合效果还可以，但不是很好！这是一个很好的基线起点，让我们能够回答曲线拟合类型的问题。我们可能不太信任超出数据范围的推断，特别是考虑到曲线对最后四个值的拟合效果不佳。例如，从50英尺处推杆的预期成功概率为：

prob_at_50 = scipy.special.expit(logit_post["a"] * 50 + logit_post["b"]).mean(dim="sample").item()
print(f"{100 * prob_at_50:.5f}%")

0.00283%

这个教训是

\[ \mathbb{E}[f(\theta)] \ne f(\mathbb{E}[\theta]). \]

这出现在这里在使用

# Right!
scipy.special.expit(logit_trace.posterior["a"] * 50 + logit_trace.posterior["b"]).mean(dim="sample")

而不是

# Wrong!
scipy.special.expit(logit_trace.posterior["a"].mean(dim="sample") * 50 + logit_trace.posterior["b"].mean(dim="sample"))

计算我们在50英尺处的期望值。

基于几何的模型

作为对这些数据进行建模的第二次尝试，既为了改进拟合度，也为了提高外推的置信度，我们考虑了这种情况的几何特性。我们假设职业高尔夫球手可以以某个方向击球，但存在一些（？）误差。具体来说，球实际行进的角度在0附近呈正态分布，我们试图学习其方差。

然后，只要角度误差足够小，球仍然能够进入杯中。这从直觉上来说很好！一个较长的推杆会允许更小的角度误差，因此成功率会比短推杆低。

我跳过了根据准确性方差和到洞口的距离计算推杆概率的推导过程，但这在几何学中是一个有趣的练习，结果是

\[ p(\text{success} | \sigma_{\text{angle}}, \text{distance}) = 2 \Phi\left( \frac{ \arcsin \left((R - r) / \text{distance}\right)}{\sigma_{\text{angle}}}\right), \]

其中 \(\Phi\) 是正态累积密度函数，\(R\) 是杯子的半径（结果是2.125英寸），\(r\) 是高尔夫球的半径（大约0.84英寸）。

为了感受这个模型，让我们来看几个手动绘制的\(\sigma_{\text{angle}}\)值。

def forward_angle_model(variances_of_shot, t):
    norm_dist = XrContinuousRV(st.norm, 0, variances_of_shot)
    return 2 * norm_dist.cdf(np.arcsin((CUP_RADIUS - BALL_RADIUS) / t)) - 1

_, ax = plt.subplots()
var_shot_ary = [0.01, 0.02, 0.05, 0.1, 0.2, 1]
var_shot_plot = xr.DataArray(var_shot_ary, coords=[("variance", var_shot_ary)])

forward_angle_model(var_shot_plot, t).plot.line(hue="variance")

plot_golf_data(golf_data, ax=ax)
ax.set_title("Model prediction for selected amounts of variance");

这看起来是一个有希望的方法！0.02弧度的方差看起来会接近正确答案。模型还预测了从0英尺开始的推杆都会进洞，这是一个不错的副作用。我们可能会考虑高尔夫球手是否对称地错过推杆。右撇子推杆手和左撇子推杆手可能会有不同的偏差，这是有可能的。

拟合模型

PyMC 已经实现了 \(\Phi\)，但它隐藏得比较深（pm.distributions.dist_math.normal_lcdf），无论如何，我们自己实现它也是值得的，使用与 误差函数 相关的恒等式。

def phi(x):
    """Calculates the standard normal cumulative distribution function."""
    return 0.5 + 0.5 * pt.erf(x / pt.sqrt(2.0))


with pm.Model() as angle_model:
    distance_ = pm.MutableData("distance", golf_data["distance"], dims="obs_id")
    tries_ = pm.MutableData("tries", golf_data["tries"], dims="obs_id")
    successes_ = pm.MutableData("successes", golf_data["successes"], dims="obs_id")

    variance_of_shot = pm.HalfNormal("variance_of_shot")
    p_goes_in = pm.Deterministic(
        "p_goes_in",
        2 * phi(pt.arcsin((CUP_RADIUS - BALL_RADIUS) / distance_) / variance_of_shot) - 1,
        dims="obs_id",
    )
    success = pm.Binomial("success", n=tries_, p=p_goes_in, observed=successes_, dims="obs_id")


pm.model_to_graphviz(angle_model)

先验预测检查

我们常常希望从先验分布中采样，特别是当我们对观测结果有一些想法，但对先验参数的直觉不多时。这里我们有一个基于角度的模型，但如果给出了角度的方差，可能不太直观地理解它如何影响射击的准确性。让我们检查一下！

有时，自定义的可视化或仪表板对于先验预测检查非常有用。在这里，我们绘制了从20英尺外推杆的先验分布。

with angle_model:
    angle_trace = pm.sample_prior_predictive(500)

angle_prior = angle_trace.prior.squeeze()

angle_of_shot = XrContinuousRV(st.norm, 0, angle_prior["variance_of_shot"]).rvs(
    random_state=RANDOM_SEED
)  # radians
distance = 20  # feet

end_positions = xr.Dataset(
    {"endx": distance * np.cos(angle_of_shot), "endy": distance * np.sin(angle_of_shot)}
)

fig, ax = plt.subplots()
for draw in end_positions["draw"]:
    end = end_positions.sel(draw=draw)
    ax.plot([0, end["endx"]], [0, end["endy"]], "k-o", lw=1, mfc="w", alpha=0.5)
ax.plot(0, 0, "go", label="Start", mfc="g", ms=20)
ax.plot(distance, 0, "ro", label="Goal", mfc="r", ms=20)

ax.set_title(f"Prior distribution of putts from {distance}ft away")

ax.legend();

Sampling: [success, variance_of_shot]

这有点有趣！最明显的是，可能不应该这么常见地把球向后推。这也让我们担心我们正在使用正态分布来建模角度。von Mises分布可能在这里是合适的。此外，高尔夫球手需要站在某个地方，所以也许为von Mises分布添加一些边界是合适的。不过，我们会发现这个模型从数据中学得相当好，这些添加并不是必要的。

with angle_model:
    angle_trace.extend(pm.sample(1000, tune=1000, target_accept=0.85))

angle_post = az.extract(angle_trace)

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

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

ax = plot_golf_data(golf_data)

angle_post["expit"] = forward_angle_model(angle_post["variance_of_shot"], t)

ax.plot(
    t,
    angle_post["expit"][:, ::100],
    lw=1,
    color="C1",
    alpha=0.1,
)

ax.plot(
    t,
    angle_post["expit"].mean(dim="sample"),
    label="Geometry-based model",
)

ax.plot(
    t,
    logit_post["expit"].mean(dim="sample"),
    label="Logit-binomial model",
)
ax.set_title("Comparing the fit of geometry-based and logit-binomial model")
ax.legend();

这个新模型似乎更符合实际情况，通过模拟情况的几何形状，我们可能对外推数据更有信心。这个模型表明，50英尺的推杆进球的几率要高得多：

angle_prob_at_50 = forward_angle_model(angle_post["variance_of_shot"], np.array([50]))
print(f"{100 * angle_prob_at_50.mean().item():.2f}% vs {100 * prob_at_50:.5f}%")

6.40% vs 0.00283%

我们也可以重新创建我们之前的预测图，这让我们对先验没有导致后验分布中出现不合理情况感到一些信心：角度的方差非常小！

angle_of_shot = XrContinuousRV(st.norm, 0, angle_post["variance_of_shot"]).rvs(
    random_state=RANDOM_SEED
)  # radians
distance = 20  # feet

end_positions = xr.Dataset(
    {"endx": distance * np.cos(angle_of_shot.data), "endy": distance * np.sin(angle_of_shot.data)}
)

fig, ax = plt.subplots()
for x, y in zip(end_positions.endx, end_positions.endy):
    ax.plot([0, x], [0, y], "k-o", lw=1, mfc="w", alpha=0.5)
ax.plot(0, 0, "go", label="Start", mfc="g", ms=20)
ax.plot(distance, 0, "ro", label="Goal", mfc="r", ms=20)

ax.set_title(f"Prior distribution of putts from {distance}ft away")
ax.set_xlim(-21, 21)
ax.set_ylim(-21, 21)
ax.legend();

新数据！

马克·布罗迪使用了新的推杆总结数据来拟合一个新的模型。我们将使用这些新数据来改进我们的模型：

#  golf putting data from Broadie (2018)
new_golf_data = """distance tries successes
0.28 45198 45183
0.97 183020 182899
1.93 169503 168594
2.92 113094 108953
3.93 73855 64740
4.94 53659 41106
5.94 42991 28205
6.95 37050 21334
7.95 33275 16615
8.95 30836 13503
9.95 28637 11060
10.95 26239 9032
11.95 24636 7687
12.95 22876 6432
14.43 41267 9813
16.43 35712 7196
18.44 31573 5290
20.44 28280 4086
21.95 13238 1642
24.39 46570 4767
28.40 38422 2980
32.39 31641 1996
36.39 25604 1327
40.37 20366 834
44.38 15977 559
48.37 11770 311
52.36 8708 231
57.25 8878 204
63.23 5492 103
69.18 3087 35
75.19 1742 24"""

new_golf_data = pd.read_csv(io.StringIO(new_golf_data), sep=" ")

ax = plot_golf_data(new_golf_data)
plot_golf_data(golf_data, ax=ax, color="C1")

t_ary = np.linspace(CUP_RADIUS - BALL_RADIUS, new_golf_data.distance.max(), 200)
t = xr.DataArray(t_ary, coords=[("distance", t_ary)])

ax.plot(
    t, forward_angle_model(angle_trace.posterior["variance_of_shot"], t).mean(("chain", "draw"))
)
ax.set_title("Comparing the new data set to the old data set, and\nconsidering the old model fit");

这个新数据集代表了大约是旧数据200倍的推杆尝试次数，并且包括了长达75英尺的推杆，而旧数据集只包括了20英尺的推杆。此外，新数据似乎代表了与旧数据不同的群体：尽管两者有不同的区间，但新数据表明大多数数据的成功率更高。这可能是由于数据收集方法的不同，或者是高尔夫球手在这几年间的进步。

在新数据上拟合模型

由于我们认为这可能是两个不同的群体，最简单的解决方案是重新拟合我们的模型。这比之前的情况更糟：出现了分歧，运行时间也更长。这可能表明模型存在问题：Andrew Gelman称之为“统计计算的民间定理”。

with angle_model:
    pm.set_data(
        {
            "distance": new_golf_data["distance"],
            "tries": new_golf_data["tries"],
            "successes": new_golf_data["successes"],
        }
    )
    new_angle_trace = pm.sample(1000, tune=1500)

100.00% [10000/10000 00:14<00:00 Sampling 4 chains, 2,003 divergences]

Sampling 4 chains for 1_500 tune and 1_000 draw iterations (6_000 + 4_000 draws total) took 15 seconds.

注意

正如您将在下面的图中看到的，这个模型对新数据的拟合效果非常差。在这种情况下，所有的分歧和收敛警告除了使用一个能够真正解释数据的模型外，没有其他解决方案。

ax = plot_golf_data(new_golf_data)
plot_golf_data(golf_data, ax=ax, color="C1")

new_angle_post = az.extract(new_angle_trace)
ax.plot(
    t,
    forward_angle_model(angle_post["variance_of_shot"], t).mean(dim="sample"),
    label="Trained on original data",
)
ax.plot(
    t,
    forward_angle_model(new_angle_post["variance_of_shot"], t).mean(dim="sample"),
    label="Trained on new data",
)
ax.set_title("Retraining the model on new data")
ax.legend();

一个包含孔洞距离的模型

我们可能会假设，除了将球击向正确的方向外，高尔夫球手还需要将球击出适当的距离。具体来说，我们假设：

1. 如果推杆短了或超过洞口3英尺以上，它就不会进洞。

2. 高尔夫球手的目标是比洞口多一英尺

3. 球的行进距离，\(u\)，是根据

   \[ u \sim \mathcal{N}\left(1 + \text{距离}, \sigma_{\text{距离}} (1 + \text{距离})\right), \]
   其中我们将学习\(\sigma_{\text{距离}}\)。

再次强调，这是一个几何和代数问题，用于计算从任意给定距离进球的概率：

\[ P(\text{好距离}) = P(\text{距离} < u < \text{距离} + 3) \]

它使用了phi，即我们之前实现的累积正态密度函数。

OVERSHOT = 1.0
DISTANCE_TOLERANCE = 3.0


with pm.Model() as distance_angle_model:
    distance_ = pm.MutableData("distance", new_golf_data["distance"], dims="obs_id")
    tries_ = pm.MutableData("tries", new_golf_data["tries"], dims="obs_id")
    successes_ = pm.MutableData("successes", new_golf_data["successes"], dims="obs_id")

    variance_of_shot = pm.HalfNormal("variance_of_shot")
    variance_of_distance = pm.HalfNormal("variance_of_distance")
    p_good_angle = pm.Deterministic(
        "p_good_angle",
        2 * phi(pt.arcsin((CUP_RADIUS - BALL_RADIUS) / distance_) / variance_of_shot) - 1,
        dims="obs_id",
    )
    p_good_distance = pm.Deterministic(
        "p_good_distance",
        phi((DISTANCE_TOLERANCE - OVERSHOT) / ((distance_ + OVERSHOT) * variance_of_distance))
        - phi(-OVERSHOT / ((distance_ + OVERSHOT) * variance_of_distance)),
        dims="obs_id",
    )

    success = pm.Binomial(
        "success", n=tries_, p=p_good_angle * p_good_distance, observed=successes_, dims="obs_id"
    )


pm.model_to_graphviz(distance_angle_model)

这个模型仍然只有2个维度需要拟合。我们可能会考虑检查OVERSHOT和DISTANCE_TOLERANCE。检查第一个可能需要给当地的高尔夫球场打电话，而第二个可能需要去果岭并进行一些实验。我们还可能考虑添加一些显式的相关性：控制角度较差可能对应于控制距离较差，或者长推杆会导致角度上的更多变化，这些都是合理的。

拟合距离角度模型

with distance_angle_model:
    distance_angle_trace = pm.sample(1000, tune=1000, target_accept=0.85)

100.00% [8000/8000 00:03<00:00 Sampling 4 chains, 0 divergences]

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

def forward_distance_angle_model(variance_of_shot, variance_of_distance, t):
    rv = XrContinuousRV(st.norm, 0, 1)
    angle_prob = 2 * rv.cdf(np.arcsin((CUP_RADIUS - BALL_RADIUS) / t) / variance_of_shot) - 1

    distance_prob_one = rv.cdf(
        (DISTANCE_TOLERANCE - OVERSHOT) / ((t + OVERSHOT) * variance_of_distance)
    )
    distance_prob_two = rv.cdf(-OVERSHOT / ((t + OVERSHOT) * variance_of_distance))
    distance_prob = distance_prob_one - distance_prob_two

    return angle_prob * distance_prob


ax = plot_golf_data(new_golf_data)

distance_angle_post = az.extract(distance_angle_trace)

ax.plot(
    t,
    forward_angle_model(new_angle_post["variance_of_shot"], t).mean(dim="sample"),
    label="Just angle",
)
ax.plot(
    t,
    forward_distance_angle_model(
        distance_angle_post["variance_of_shot"],
        distance_angle_post["variance_of_distance"],
        t,
    ).mean(dim="sample"),
    label="Distance and angle",
)

ax.set_title("Comparing fits of models on new data")
ax.legend();

这个新模型看起来更好，并且与旧模型相比，在采样问题更少的情况下更快地拟合。在10到40英尺之间存在一些不匹配，但总体上看起来不错。我们可以通过获取后验预测样本并查看残差来得出相同的结论。在这里，我们看到拟合是由前4个区间驱动的，这些区间包含了大约40%的数据。

with distance_angle_model:
    pm.sample_posterior_predictive(distance_angle_trace, extend_inferencedata=True)

const_data = distance_angle_trace.constant_data
pp = distance_angle_trace.posterior_predictive
residuals = 100 * ((const_data["successes"] - pp["success"]) / const_data["tries"]).mean(
    ("chain", "draw")
)

fig, ax = plt.subplots()

ax.plot(new_golf_data.distance, residuals, "o-")
ax.axhline(y=0, linestyle="dashed", linewidth=1)
ax.set_xlabel("Distance from hole")
ax.set_ylabel("Absolute error in expected\npercent of success")

ax.set_title("Residuals of new model");

100.00% [4000/4000 00:00<00:00]

一个新模型

在此处停止是合理的，但如果我们想在各处改进拟合效果，我们可能希望选择一个不同于Binomial的似然函数，后者非常关注那些有大量观测点的数据。我们可以做的一件事是为每个数据点添加一些独立的额外误差。我们可以通过几种方式来实现这一点：

1. 二项分布通常由\(n\)（观察次数）和\(p\)（单次成功的概率）参数化。我们可以改为由均值（\(np\)）和方差（\(np(1-p)\)）参数化，并在似然函数中添加与\(n\)无关的误差。

2. 使用 BetaBinomial 分布，尽管那里的误差仍然会（大致）与观察数量成正比

3. 用正态分布近似二项分布的成功概率。这实际上等同于第一种方法，但不需要自定义分布。注意，我们将使用\(p\)作为均值，并使用\(p(1-p) / n\)作为方差。一旦我们加入一些离散度\(\epsilon\)，方差变为\(p(1-p)/n + \epsilon\)。

我们遵循方法3，如在Stan案例研究中所示，并将1留作练习。

with pm.Model() as disp_distance_angle_model:
    distance_ = pm.MutableData("distance", new_golf_data["distance"], dims="obs_id")
    tries_ = pm.MutableData("tries", new_golf_data["tries"], dims="obs_id")
    successes_ = pm.MutableData("successes", new_golf_data["successes"], dims="obs_id")
    obs_prop_ = pm.MutableData(
        "obs_prop", new_golf_data["successes"] / new_golf_data["tries"], dims="obs_id"
    )

    variance_of_shot = pm.HalfNormal("variance_of_shot")
    variance_of_distance = pm.HalfNormal("variance_of_distance")
    dispersion = pm.HalfNormal("dispersion")

    p_good_angle = pm.Deterministic(
        "p_good_angle",
        2 * phi(pt.arcsin((CUP_RADIUS - BALL_RADIUS) / distance_) / variance_of_shot) - 1,
        dims="obs_id",
    )
    p_good_distance = pm.Deterministic(
        "p_good_distance",
        phi((DISTANCE_TOLERANCE - OVERSHOT) / ((distance_ + OVERSHOT) * variance_of_distance))
        - phi(-OVERSHOT / ((distance_ + OVERSHOT) * variance_of_distance)),
        dims="obs_id",
    )

    p = p_good_angle * p_good_distance
    p_success = pm.Normal(
        "p_success",
        mu=p,
        sigma=pt.sqrt(((p * (1 - p)) / tries_) + dispersion**2),
        observed=obs_prop_,  # successes_ / tries_
        dims="obs_id",
    )


pm.model_to_graphviz(disp_distance_angle_model)

with disp_distance_angle_model:
    disp_distance_angle_trace = pm.sample(1000, tune=1000)
    pm.sample_posterior_predictive(disp_distance_angle_trace, extend_inferencedata=True)

100.00% [4000/4000 00:00<00:00]

ax = plot_golf_data(new_golf_data, ax=None)

disp_distance_angle_post = az.extract(disp_distance_angle_trace)

ax.plot(
    t,
    forward_distance_angle_model(
        distance_angle_post["variance_of_shot"],
        distance_angle_post["variance_of_distance"],
        t,
    ).mean(dim="sample"),
    label="Distance and angle",
)
ax.plot(
    t,
    forward_distance_angle_model(
        disp_distance_angle_post["variance_of_shot"],
        disp_distance_angle_post["variance_of_distance"],
        t,
    ).mean(dim="sample"),
    label="Dispersed model",
)
ax.set_title("Comparing dispersed model with binomial distance/angle model")
ax.legend();

这个新模型在10到30英尺之间表现更好，正如我们通过残差图也能看到的那样——请注意，这个模型在非常短的推杆上表现略差：

const_data = distance_angle_trace.constant_data
old_pp = az.extract(distance_angle_trace, group="posterior_predictive")
old_residuals = 100 * ((const_data["successes"] - old_pp["success"]) / const_data["tries"]).mean(
    dim="sample"
)

pp = az.extract(disp_distance_angle_trace, group="posterior_predictive")
residuals = 100 * (const_data["successes"] / const_data["tries"] - pp["p_success"]).mean(
    dim="sample"
)

fig, ax = plt.subplots()

ax.plot(new_golf_data.distance, residuals, label="Dispersed model")
ax.plot(new_golf_data.distance, old_residuals, label="Distance and angle model")
ax.legend()
ax.axhline(y=0, linestyle="dashed", linewidth=1)
ax.set_xlabel("Distance from hole")
ax.set_ylabel("Absolute error in expected\npercent of success")
ax.set_title("Residuals of dispersed model vs distance/angle model");

超越预测

我们希望使用贝叶斯分析，因为我们关心量化参数中的不确定性。我们有一个漂亮的几何模型，它不仅为我们提供预测，还为我们提供了参数的后验分布。我们可以利用这一点来推断如果不在洞中，我们的推杆可能落在哪里！

首先，我们可以尝试可视化一位职业高尔夫球手的20,000次推杆可能的样子。我们：

1. 将试验次数设置为5

2. 对于每个联合后验样本 variance_of_shot 和 variance_of_distance， 从正态分布中抽取一个角度和一个距离5次。

3. 绘制点，除非它会落入洞中

def simulate_from_distance(trace, distance_to_hole, trials=5):
    variance_of_shot = trace.posterior["variance_of_shot"]
    variance_of_distance = trace.posterior["variance_of_distance"]

    theta = XrContinuousRV(st.norm, 0, variance_of_shot).rvs(size=trials, dims="trials")
    distance = XrContinuousRV(
        st.norm, distance_to_hole + OVERSHOT, (distance_to_hole + OVERSHOT) * variance_of_distance
    ).rvs(size=trials, dims="trials")

    final_position = xr.concat(
        (distance * np.cos(theta), distance * np.sin(theta)), dim="axis"
    ).assign_coords(axis=["x", "y"])

    made_it = np.abs(theta) < np.arcsin((CUP_RADIUS - BALL_RADIUS) / distance_to_hole)
    made_it = (
        made_it
        * (final_position.sel(axis="x") > distance_to_hole)
        * (final_position.sel(axis="x") < distance_to_hole + DISTANCE_TOLERANCE)
    )

    dims = [dim for dim in final_position.dims if dim != "axis"]
    final_position = final_position.where(~made_it).stack(idx=dims).dropna(dim="idx")
    total_simulations = made_it.size

    _, ax = plt.subplots()

    ax.plot(0, 0, "k.", lw=1, mfc="black", ms=250 / distance_to_hole)
    ax.plot(*final_position, ".", alpha=0.1, mfc="r", ms=250 / distance_to_hole, mew=0.5)
    ax.plot(distance_to_hole, 0, "ko", lw=1, mfc="black", ms=350 / distance_to_hole)

    ax.set_facecolor("#e6ffdb")
    ax.set_title(
        f"Final position of {total_simulations:,d} putts from {distance_to_hole}ft.\n"
        f"({100 * made_it.mean().item():.1f}% made)"
    )
    return ax


simulate_from_distance(distance_angle_trace, distance_to_hole=50);

simulate_from_distance(distance_angle_trace, distance_to_hole=7);

然后，我们可以利用这一点来计算球员从给定距离可能需要进行的推杆次数。这可以影响战略决策，例如尝试以更少的击球次数到达果岭，这可能导致更长的第一次推杆，与更保守的方法相比。我们通过模拟推杆直到它们全部进洞来实现这一点。

请注意，这又是我们可能会通过实验验证的内容。特别是，在互联网上进行的一项非常不科学的搜索发现，专业人士在20-25英尺的距离内只有大约3%的时间会进行3推。我们的模型预测在22.5英尺的距离内进行3次或更多推杆的概率为2.8%，这看起来好得令人怀疑。

def expected_num_putts(trace, distance_to_hole, trials=100_000):
    distance_to_hole = distance_to_hole * np.ones(trials)

    combined_trace = az.extract(trace)

    n_samples = combined_trace.dims["sample"]

    idxs = np.random.randint(0, n_samples, trials)
    variance_of_shot = combined_trace["variance_of_shot"].isel(sample=idxs)
    variance_of_distance = combined_trace["variance_of_distance"].isel(sample=idxs)
    n_shots = []
    while distance_to_hole.size > 0:
        theta = np.random.normal(0, variance_of_shot)
        distance = np.random.normal(
            distance_to_hole + OVERSHOT, (distance_to_hole + OVERSHOT) * variance_of_distance
        )

        final_position = np.array([distance * np.cos(theta), distance * np.sin(theta)])

        made_it = np.abs(theta) < np.arcsin(
            (CUP_RADIUS - BALL_RADIUS) / distance_to_hole.clip(min=CUP_RADIUS - BALL_RADIUS)
        )
        made_it = (
            made_it
            * (final_position[0] > distance_to_hole)
            * (final_position[0] < distance_to_hole + DISTANCE_TOLERANCE)
        )

        distance_to_hole = np.sqrt(
            (final_position[0] - distance_to_hole) ** 2 + final_position[1] ** 2
        )[~made_it].copy()
        variance_of_shot = variance_of_shot[~made_it]
        variance_of_distance = variance_of_distance[~made_it]
        n_shots.append(made_it.sum())
    return np.array(n_shots) / trials

distances = (10, 20, 40, 80)
fig, axes = plt.subplots(nrows=2, ncols=2, sharex=True, sharey=True, figsize=(10, 10))

for distance, ax in zip(distances, axes.ravel()):
    made = 100 * expected_num_putts(disp_distance_angle_trace, distance)
    x = np.arange(1, 1 + len(made), dtype=int)
    ax.vlines(np.arange(1, 1 + len(made)), 0, made, linewidths=50)
    ax.set_title(f"{distance} feet")
    ax.set_ylabel("Percent of attempts")
    ax.set_xlabel("Number of putts")
ax.set_xticks(x)
ax.set_ylim(0, 100)
ax.set_xlim(0, 5.6)
fig.suptitle("Simulated number of putts from\na few distances");

作者

• 由Colin Carroll改编自Stan文档中的[高尔夫推杆的模型构建与扩展]案例研究（pymc#3666）

• 由 Marco Gorelli 更新（pymc-examples#39）

• 由 Oriol Abril-Pla 更新以使用 PyMC v4 和 xarray-einstats

• 由Benjamin T. Vincent更新，以在2023年2月使用az.extract（pymc-examples#522）

参考资料

[1]

唐纳德·A·贝里。《统计学：贝叶斯视角》。Duxbury出版社，1996年。

水印

%load_ext watermark
%watermark -n -u -v -iv -w -p aeppl,xarray_einstats

Last updated: Sun Feb 05 2023

Python implementation: CPython
Python version       : 3.10.9
IPython version      : 8.9.0

aeppl          : not installed
xarray_einstats: 0.5.1

xarray    : 2023.1.0
arviz     : 0.14.0
numpy     : 1.24.1
matplotlib: 3.6.3
pymc      : 5.0.1
scipy     : 1.10.0
pandas    : 1.5.3
pytensor  : 2.8.11

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"
}

渲染后可能看起来像：