C01 · Golf putting: physics beats curve fitting — Bayesian Modelling in Practice

The brief

A golf analytics company wants a model for the probability that a professional sinks a putt as a function of distance. It will feed a "strokes gained" product, so it has to extrapolate sensibly and its parameters should mean something to a coach.

You get two datasets: a small classic one (19 distances, ~5,600 putts) and a modern one with 30x more putts. A model that looks great on the first will break on the second. Your job is to find out why and fix it.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task2", sigma_angle_deg=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
from scipy import stats

from pymc_challenges import Hints, data

RANDOM_SEED = 1996
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C01")
h.tasks()

C01 - Golf putting: when physics beats curve fitting

task1 - Logistic regression baseline (3 hints, 2 check(s))
task2 - A mechanistic angle model (3 hints, 1 check(s))
task3 - LOO comparison (3 hints)
task4 - New data breaks the model (3 hints)
task5 - Distance control and a sampler that will not mix (3 hints, 2 check(s))
task6 - Answer the business question (3 hints, 1 check(s))

The data

x is the distance to the hole in feet, n the number of putts attempted from that distance and y how many went in.

data.describe("golf")
golf = data.load("golf")
golf.head()

golf: Putts attempted and made by distance (feet), 19 distances.
  source: Berry (1996), professional golf putting; used in Gelman & Nolan (2002).
  url:    https://raw.githubusercontent.com/stan-dev/example-models/master/knitr/golf/golf_data.txt

	x	n	y
0	2	1443	1346
1	3	694	577
2	4	455	337
3	5	353	208
4	6	272	149

Two physical constants you will need later (both in feet, like x):

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

Task 0 · Look at the data

Plot the empirical success rate against distance with an uncertainty interval for each point. Which distances are measured precisely, and which are not?

Your turn Work this task in the notebook.

Task 1 · The default: logistic regression

Fit the model every analyst reaches for first:

$$y_i \sim \text{Binomial}(n_i, p_i), \qquad \text{logit}(p_i) = a + b\,x_i$$

Deliver

Priors you can defend. Show with a prior predictive plot that they imply plausible curves (no 50% success at 0 feet, no cliff edges).
A fitted model with clean diagnostics (r_hat, ESS, divergences).
A plot of the posterior success curve with a credible band on top of the data. Where does the fit fail, and in which direction?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

You have counts of successes out of n attempts at each distance, not individual putts. Which likelihood is that? And before sampling anything: if you draw a and b from your priors and plot the implied curves over 0-20 feet, do they look like things a golfer could do?

2 Approach

Use pm.Binomial("y", n=..., logit_p=a + b * x, observed=...) - the logit_p argument saves you from writing the inverse-logit yourself and is numerically safer.

For the prior predictive check, pm.sample_prior_predictive(200) returns a DataTree with a prior group. az.extract(idata, group="prior", var_names=["a", "b"]) gives you flat draws you can broadcast against a grid of distances in plain NumPy.

On the prior scale: distance spans ~20 feet, so a slope of 0.5 per foot is already a 10-unit swing in log-odds. Priors much wider than that produce step functions.

3 Skeleton

with pm.Model() as logit_model:
    a = pm.Normal("a", 0, ___)
    b = pm.Normal("b", 0, ___)
    pm.Binomial("y", n=golf.n.values, logit_p=___, observed=golf.y.values)
    prior = pm.sample_prior_predictive(200, random_seed=RANDOM_SEED)
    idata = pm.sample(random_seed=RANDOM_SEED)

post = az.extract(idata, var_names=["a", "b"])          # dims: (sample,)
logit = post["a"].values[:, None] + post["b"].values[:, None] * x_grid
curves = 1 / (1 + np.exp(-logit))                       # (draws, grid)
lo, hi = np.quantile(curves, [0.03, 0.97], axis=0)

Self-checks in the notebook (h.check("task1", …)):

a — posterior mean of the intercept, roughly 2.2 on the logit scale
b — posterior mean of the slope, roughly -0.26 per foot

Task 2 · Think like a golfer

A putt goes in if the ball's centre passes over the hole. Suppose the golfer aims at the centre but the launch angle is off by a random error $\theta \sim \text{Normal}(0, \sigma)$.

Deliver

Derive $p(x)$ from the geometry (constants BALL_RADIUS, HOLE_RADIUS are above) and write it down in a markdown cell. It has one free parameter.
Implement it in PyMC and fit it. Report $\sigma$ in degrees.
Overlay its posterior curve on the logistic fit. One parameter versus two - which wins?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Draw it. The ball starts x feet from the hole. For its centre to pass over the hole, how far off-line can it be when it gets there - in terms of the hole radius R and ball radius r? That gives you a right-angled triangle, and the triangle gives you a threshold angle.

2 Approach

The threshold angle is theta_0 = arcsin((R - r) / x). If the angular error is Normal(0, sigma), the putt drops when |theta| < theta_0, so

p = 2 * Phi(theta_0 / sigma) - 1

The standard normal CDF Phi is pm.math.invprobit in PyMC; pm.math.arcsin exists too. sigma is in radians and must be positive - think about how many degrees a professional could plausibly be off by when choosing the prior scale. Wrap the conversion to degrees in a pm.Deterministic so it shows up in az.summary.

3 Skeleton

with pm.Model() as angle_model:
    sigma_angle = pm.HalfNormal("sigma_angle", ___)          # radians
    pm.Deterministic("sigma_angle_deg", sigma_angle * 180 / np.pi)
    theta_0 = pm.math.arcsin((HOLE_RADIUS - BALL_RADIUS) / golf.x.values)
    p = 2 * pm.math.invprobit(___) - 1
    pm.Binomial("y", n=golf.n.values, p=p, observed=golf.y.values)
    idata = pm.sample(random_seed=RANDOM_SEED)

For plotting, re-implement p(x) with NumPy and scipy.stats.norm.cdf and evaluate it for every posterior draw of sigma_angle on a grid.

Self-checks in the notebook (h.check("task2", …)):

sigma_angle_deg — posterior mean of the angular error in DEGREES, about 1.5 (did you convert from radians?)

Task 3 · Make the comparison quantitative

Compare the two models with PSIS-LOO.

Deliver the az.compare table and two sentences: what does it say, and how far do you trust it here? (Look at the warnings - there are only 19 "observations".)

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

LOO needs the pointwise log-likelihood of each observation under each posterior draw. pm.sample does not store it by default. How do you get it after the fact?

2 Approach

pm.compute_log_likelihood(idata) inside the model context adds a log_likelihood group in place. Then az.compare({"name": idata, ...}) returns a DataFrame ranked by elpd.

Read the warnings: Pareto-k above 0.7 means the importance-sampling approximation behind PSIS-LOO is unreliable for that observation. Ask yourself what "leaving one observation out" means when one observation is 1,400 putts.

3 Skeleton

with logit_model:
    pm.___(logit_idata)          # adds the log_likelihood group in place
with angle_model:
    pm.___(angle_idata)
az.compare({"logistic": ___, "angle": ___})
az.loo(angle_idata)     # prints the Pareto-k table for one model

Task 4 · New data breaks the model

Mark Broadie collected far more putts, out to 75 feet.

Deliver

Fit your angle model to golf_new, unchanged. Plot fit against data, and plot the residuals (observed rate − posterior mean). What is the pattern?
A one-paragraph hypothesis: what real-world mechanism is the model missing?

data.describe("golf_new")
golf_new = data.load("golf_new")
golf_new.tail()

golf_new: Putts attempted and made by distance (feet), 31 distances, far more putts.
  source: Mark Broadie's modern PGA putting data, via Gelman's golf case study.
  url:    https://raw.githubusercontent.com/stan-dev/example-models/master/knitr/golf/golf_data_new.txt

	x	n	y
26	52.36	8708	231
27	57.25	8878	204
28	63.23	5492	103
29	69.18	3087	35
30	75.19	1742	24

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Refit without changing anything, then look at where the misfit is rather than how big it is. Is the sign of the residual random, or does it depend on distance? Note also what happened to the estimate of sigma compared with Task 2.

2 Approach

Residual = observed rate - posterior mean of p at each observed distance. You will see the model is too optimistic for long putts. Angle is not the only way to miss: think about what else has to be right for a 50-foot putt to drop.

3 Skeleton

sigma = az.extract(idata_new, var_names="sigma_angle").values        # (draws,)
p_hat = p_angle_numpy(golf_new.x.values, sigma[:, None]).mean(axis=___)   # average over draws
plt.plot(golf_new.x, ___ - p_hat, "o-")                                   # observed rate minus fitted
plt.axhline(0)

The missing mechanism is distance control: hit it too softly and it stops short, too hard and it hops over the hole.

Task 5 · Add distance control (and survive the sampler)

Domain knowledge from a putting coach:

Golfers aim to hit the ball about 1 foot past the hole.
The putt can drop if the ball would have travelled between the hole and 3 feet past it.
The error in how far the ball travels is proportional to the intended distance.

Deliver

Turn that into $p_\text{distance}(x)$ with one new parameter, and combine it with the angle component (assume independence).
Fit it with a Binomial likelihood. Check the diagnostics carefully. Something will be wrong - explain why this posterior is hard, in terms of the data.
Fix the model (not the sampler settings!) and refit. Report all parameters and show the final fit and residuals.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Part 1: the golfer intends the ball to travel x + 1 feet. If the error is proportional to that intended distance, what is the distribution of the distance actually travelled? The putt drops if that lands in an interval - which one?

Part 2: when the chains disagree, do not reach for target_accept. Look at the n column. How precisely does a Binomial with n = 50,000 claim to know each success rate? Can any smooth two-parameter curve pass through 31 points that precisely?

2 Approach

Distance travelled is (x + 1) * (1 + eps) with eps ~ Normal(0, sigma_distance), and it must land in [x, x + 3]:

p_distance = Phi(2 / ((x + 1) * sigma_d)) - Phi(-1 / ((x + 1) * sigma_d))

and p = p_angle * p_distance.

The Binomial version is multimodal because the model is slightly wrong and the data are abundant enough to punish that brutally. The honest fix is to add a term for model error: let the observed rate y / n be Normal around p with variance equal to the binomial sampling variance p (1 - p) / n plus an extra sigma_y ** 2.

3 Skeleton

with pm.Model() as final_model:
    sigma_angle = pm.HalfNormal("sigma_angle", 0.1)
    sigma_distance = pm.HalfNormal("sigma_distance", ___)
    sigma_y = pm.HalfNormal("sigma_y", ___)          # on the probability scale: small!
    scale = (x + 1) * sigma_distance
    p_dist = pm.math.invprobit(2 / scale) - pm.math.invprobit(___)
    p = p_angle * p_dist
    pm.Normal("rate", mu=p,
              sigma=pm.math.sqrt(p * (1 - p) / n + sigma_y**2),
              observed=y / n)

Self-checks in the notebook (h.check("task5", …)):

sigma_angle_deg — posterior mean of angular error in degrees in the FINAL model, about 1.0
sigma_distance — posterior mean of the proportional distance error in the FINAL model, about 0.08

Task 6 · Answer the business question

The product team asks: "From what distance is a tour pro a coin flip?"

Deliver the posterior distribution of $x_{50}$, the distance at which $p(x) = 0.5$, with a point estimate and a 94% interval. Then the follow-up: a coaching intervention claims to cut distance error $\sigma_d$ by a fifth. How many extra putts per 100 would that sink from 30 feet?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

x50 is a deterministic function of the parameters, so every posterior draw of the parameters gives you one draw of x50. There is no closed form - but you do not need one.

2 Approach

Evaluate p(x) on a fine grid of distances for every posterior draw (a (draws, grid) array via broadcasting), then for each row find the first grid point where p < 0.5. np.argmax(condition, axis=1) returns the index of the first True per row.

The what-if is the same trick: compute p(30) per draw with sigma_distance and again with 0.8 * sigma_distance, and summarise the difference.

3 Skeleton

post = az.extract(final_idata, var_names=["sigma_angle", "sigma_distance"])
sa = post["sigma_angle"].values[:, None]
sd = post["sigma_distance"].values[:, None]
grid = np.linspace(1, 40, 4000)
p = ___                                                       # (draws, grid), by broadcasting
x50 = grid[np.argmax(___, axis=1)]                            # first grid point where ... per draw
np.quantile(x50, [0.03, 0.97])

Self-checks in the notebook (h.check("task6", …)):

x50 — posterior mean distance (feet) at which p = 0.5 under the final model, about 8

Going further

The final residuals still wiggle systematically between 3 and 12 feet. Propose and test a mechanism.
Refit the final model with the exact Binomial likelihood plus a Beta-Binomial or a logit-normal random effect per distance instead of the Normal approximation. Do the conclusions change?
Fit both datasets jointly with shared physics but dataset-specific skill parameters. Did pros improve between the 1990s and the 2010s?

h.progress()

C01: 0/18 hints used, 6/6 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Look at the data

#Task 1 · The default: logistic regression

#Task 2 · Think like a golfer

#Task 3 · Make the comparison quantitative

#Task 4 · New data breaks the model

#Task 5 · Add distance control (and survive the sampler)

#Task 6 · Answer the business question

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Look at the data

Task 1 · The default: logistic regression

Task 2 · Think like a golfer

Task 3 · Make the comparison quantitative

Task 4 · New data breaks the model

Task 5 · Add distance control (and survive the sampler)

Task 6 · Answer the business question

Going further