In [1]:
import io

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy
import scipy.stats as st

print('Running on PyMC3 v{}'.format(pm.__version__))
Running on PyMC3 v3.9.0
In [2]:
plt.rcParams.update({
    "axes.prop_cycle": plt.cycler("color", ['#000000', '#1b6989', '#e69f00', '#009e73', '#f0e442', '#50b4e9', '#d55e00', '#cc79a7']),
    "figure.figsize": [12.0, 5.0],
    "font.serif": ['Palatino',
                   'Palatino Linotype',
                   'Palatino LT STD',
                   'Book Antiqua',
                   'Georgia',
                   'DejaVu Serif'],
    'font.family': 'serif',
    'figure.facecolor': '#fffff8',
    'figure.constrained_layout.use': True,
    'font.size': 14.0,
    'hist.bins': 'auto',
    'lines.linewidth': 3.0,
    'lines.markeredgewidth': 2.0,
    'lines.markerfacecolor': 'none',
    'lines.markersize': 8.0, 
})

Model building and expansion for golf putting

This uses and closely follows the case study from Andrew Gelman, written in Stan. There are some new visualizations and I steered away from using improper priors, but much credit to him and to the Stan group for the wonderful case study and software.

We use a data set from "Statistics: A Bayesian Perspective", by Don Berry (1995). The dataset describes the outcome of professional golfers putting from a number of distances, and is small enough that we can just print and load it inline, instead of doing any special csv reading.

In [3]:
# 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=" ")
In [4]:
def plot_golf_data(golf_data, ax=None):
    """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)
    ax.plot(golf_data.distance, golf_data.successes / golf_data.tries, 'o', 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 (1995)");

After plotting, we see that generally golfers are less accurate from further away. Note that this data is pre-aggregated: we may be able to do more interesting work with granular putt-by-putt data. This data set appears to have been binned to the nearest foot.

We might think about doing prediction with this data: fitting a curve to this data would allow us to make reasonable guesses at intermediate distances, as well as perhaps to extrapolate to longer distances.

Logit model

First we will fit a traditional logit-binomial model. We model the number of successes directly, with

$$ 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})) $$

Here is how to write that model in PyMC3. We wrap our models in functions to avoid polluting the namespace, and to let us swap out the data later, when we will work with a newer data set.

In [5]:
def logit_model(golf_data):
    with pm.Model() as logit_binomial:
        a = pm.Normal('a')
        b = pm.Normal('b')

        success = pm.Binomial('success', 
                              n=golf_data.tries, 
                              p=pm.math.invlogit(a * golf_data.distance + b), 
                              observed=golf_data.successes)
    return logit_binomial


pm.model_to_graphviz(logit_model(golf_data))
Out[5]:
%3 cluster19 19 b b ~ Normal success success ~ Binomial b->success a a ~ Normal a->success

We have some intuition that $a$ should be negative, and also that $b$ should be positive (since when $\text{distance} = 0$, we expect to make nearly 100% of putts). We are not putting that into the model, though. We are using this as a baseline, and we may as well wait and see if we need to add stronger priors.

In [6]:
with logit_model(golf_data):
    logit_trace = pm.sample(1000, tune=1000)
    

pm.summary(logit_trace)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [b, a]
100.00% [8000/8000 00:08<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 10 seconds.
The acceptance probability does not match the target. It is 0.8835381150026076, but should be close to 0.8. Try to increase the number of tuning steps.
The acceptance probability does not match the target. It is 0.8833160780005165, but should be close to 0.8. Try to increase the number of tuning steps.
The number of effective samples is smaller than 25% for some parameters.
/dependencies/arviz/arviz/data/io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  FutureWarning,
Out[6]:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_mean ess_sd ess_bulk ess_tail r_hat
a -0.255 0.007 -0.269 -0.244 0.000 0.000 899.0 895.0 906.0 982.0 1.0
b 2.224 0.060 2.112 2.338 0.002 0.001 930.0 925.0 931.0 935.0 1.0

We see $a$ and $b$ have the signs we expected. There were no bad warnings emitted from the sampler. Looking at the summary, the number of effective samples is reasonable, and the rhat is close to 1. This is a small model, so we are not being too careful about inspecting the fit.

We plot 50 posterior draws of $p(\text{success})$ along with the expected value. Also, we draw 500 points from the posterior predictive to plot:

In [7]:
# Draw posterior predictive samples
with logit_model(golf_data):
    # hard to plot more than 500 sensibly
    logit_ppc = pm.sample_posterior_predictive(logit_trace, samples=500)

logit_ppc_success = logit_ppc['success'].T / golf_data.tries.values.reshape(-1, 1)

# Plotting
ax = plot_golf_data(golf_data)
t = np.linspace(0, golf_data.distance.max(), 200)

for idx in np.random.randint(0, len(logit_trace), 50):
    ax.plot(t, scipy.special.expit(logit_trace['a'][idx] * t + logit_trace['b'][idx]), lw=1, color='C1', alpha=0.5)
ax.plot(t, scipy.special.expit(logit_trace['a'].reshape(-1, 1) * t + logit_trace['b'].reshape(-1, 1)).mean(axis=0), color='C2')

ax.plot(golf_data.distance, logit_ppc_success, 'k.', alpha=0.01)
ax.set_title('Logit mean and posterior predictive');
/dependencies/pymc3/pymc3/sampling.py:1618: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  "samples parameter is smaller than nchains times ndraws, some draws "
100.00% [500/500 00:02<00:00]

The fit is ok, but not great! It is a good start for a baseline, and lets us answer curve-fitting type questions. We may not trust much extrapolation beyond the end of the data, especially given how the curve does not fit the last four values very well. For example, putts from 50 feet are expected to be made with probability:

In [8]:
print(f"{100 * scipy.special.expit(logit_trace['a'] * 50 + logit_trace['b']).mean():.5f}%")
0.00280%

Aside: Calculating expectations with traces

An interesting thing to note about calculating the posterior expectation -- a strength of Bayesian models is that you can get strong correlations between parameters, which often show up particularly in regression models:

In [9]:
ax = pm.plot_pair(logit_trace, figsize=(7, 5));

corr = np.correlate(logit_trace['a'], logit_trace['b'])[0]
ax.set_title(f"Correlation: {corr:,.2f}");
/dependencies/arviz/arviz/data/io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  FutureWarning,

The lesson from this is that

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

this appeared here in using

# Right!
scipy.special.expit(logit_trace['a'] * 50 + logit_trace['b']).mean()

rather than

# Wrong!
scipy.special.expit(logit_trace['a'].mean() * 50 + logit_trace['b'].mean())

to calculate our expectation at 50 feet.

Geometry-based model

As a second pass at modelling this data, both to improve fit and to increase confidence in extrapolation, we think about the geometry of the situation. We suppose professional golfers can hit the ball in a certain direction, with some small(?) error. Specifically, the angle the ball actually travels is normally distributed around 0, with some variance that we will try to learn.

Then the ball goes in whenever the error in angle is small enough that the ball still hits the cup. This is intuitively nice! A longer putt will admit a smaller error in angle, and so a lower success rate than for shorter putts.

I am skipping a derivation of the probability of making a putt given the accuracy variance and distance to the hole, but it is a fun exercise in geometry, and turns out to be

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

where $\Phi$ is the normal cumulative density function, $R$ is the radius of the cup (turns out 2.125 inches), and $r$ is the radius of the golf ball (around 0.84 inches).

To get a feeling for this model, let's look at a few manually plotted values for $\sigma_{\text{angle}}$.

In [10]:
BALL_RADIUS = (1.68 / 2) / 12
CUP_RADIUS = (4.25 / 2) / 12

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

ax = plot_golf_data(golf_data)

t = np.linspace(0, golf_data.distance.max(), 200)

for variance_of_shot in (0.01, 0.02, 0.05, 0.1, 0.2, 1):
    ax.plot(t, forward_angle_model(variance_of_shot, t), label=f"variance={variance_of_shot}")
ax.set_title("Model prediction for selected amounts of variance")
ax.legend();

This looks like a promising approach! A variance of 0.02 radians looks like it will be close to the right answer. The model also predicted that putts from 0 feet all go in, which is a nice side effect. We might think about whether a golfer misses putts symmetrically. It is plausible that a right handed putter and a left handed putter might have a different bias to their shots.

Fitting the model

PyMC3 has $\Phi$ implemented, but it is pretty hidden (pm.distributions.dist_math.normal_lcdf), and it is worthwhile to implement it ourselves anyways, using an identity with the error function.

In [11]:
import theano.tensor as tt


def Phi(x):
    """Calculates the standard normal cumulative distribution function."""
    return 0.5 + 0.5 * tt.erf(x / tt.sqrt(2.))

def angle_model(golf_data):
    with pm.Model() as angle_model:
        variance_of_shot = pm.HalfNormal('variance_of_shot')
        p_goes_in = pm.Deterministic('p_goes_in', 2 * Phi(tt.arcsin((CUP_RADIUS - BALL_RADIUS) / golf_data.distance) / variance_of_shot) - 1)
        success = pm.Binomial('success', n=golf_data.tries, p=p_goes_in, observed=golf_data.successes)
    return angle_model

pm.model_to_graphviz(angle_model(golf_data))
Out[11]:
%3 cluster19 19 variance_of_shot variance_of_shot ~ HalfNormal p_goes_in p_goes_in ~ Deterministic variance_of_shot->p_goes_in success success ~ Binomial p_goes_in->success

Prior Predictive Checks

We often wish to sample from the prior, especially if we have some idea of what the observations would look like, but not a lot of intuition for the prior parameters. We have an angle-based model here, but it might not be intuitive if the variance of the angle is given, how that effects the accuracy of a shot. Let's check!

Sometimes a custom visualization or dashboard is useful for a prior predictive check. Here, we plot our prior distribution of putts from 20 feet away.

In [12]:
with angle_model(golf_data):
    angle_prior = pm.sample_prior_predictive(500)

angle_of_shot = np.random.normal(0, angle_prior['variance_of_shot'])  # radians
distance = 20 # feet

end_positions = np.array([
    distance * np.cos(angle_of_shot),
    distance * np.sin(angle_of_shot)
])

fig, ax = plt.subplots()
for endx, endy in end_positions.T:
    ax.plot([0, endx], [0, 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();

This is a little funny! Most obviously, it should probably be not this common to putt the ball backwards. This also leads us to worry that we are using a normal distribution to model an angle. The von Mises distribution may be appropriate here. Also, the golfer needs to stand somewhere, so perhaps adding some bounds to the von Mises would be appropriate. We will find that this model learns from the data quite well, though, and these additions are not necessary.

In [13]:
with angle_model(golf_data):
    angle_trace = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [variance_of_shot]
100.00% [8000/8000 00:06<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 7 seconds.
In [14]:
ax = plot_golf_data(golf_data)

t = np.linspace(0, golf_data.distance.max(), 200)

for idx in np.random.randint(0, len(angle_trace), 50):
    ax.plot(t, forward_angle_model(angle_trace['variance_of_shot'][idx], t), lw=1, color='C1', alpha=0.1)
ax.plot(t, forward_angle_model(angle_trace['variance_of_shot'].mean(), t), label='Geometry-based model')

ax.plot(t, scipy.special.expit(logit_trace['a'].mean() * t + logit_trace['b'].mean()), label='Logit-binomial model')
ax.set_title("Comparing the fit of geometry-based and logit-binomial model")
ax.legend();

This new model appears to fit much better, and by modelling the geometry of the situation, we may have a bit more confidence in extrapolating the data. This model suggests that a 50 foot putt has much higher chance of going in:

In [15]:
print(f"{100 * forward_angle_model(angle_trace['variance_of_shot'], 50).mean():.2f}% vs {100 * scipy.special.expit(logit_trace['a'] * 50 + logit_trace['b']).mean():.5f}%")
6.41% vs 0.00280%

We can also recreate our prior predictive plot, giving us some confidence that the prior was not leading to unreasonable situations in the posterior distribution: the variance in angle is quite small!

In [16]:
angle_of_shot = np.random.normal(0, np.random.choice(angle_trace['variance_of_shot'].flatten(), 500))  # radians
distance = 20 # feet

end_positions = np.array([
    distance * np.cos(angle_of_shot),
    distance * np.sin(angle_of_shot)
])

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

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

New Data!

Mark Broadie used new summary data on putting to fit a new model. We will use this new data to refine our model:

In [17]:
#  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=" ")
In [18]:
ax = plot_golf_data(new_golf_data)
plot_golf_data(golf_data, ax=ax)
t = np.linspace(0, new_golf_data.distance.max(), 200)

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

This new data set represents ~200 times the number of putt attempts as the old data, and includes putts up to 75ft, compared to 20ft for the old data set. It also seems that the new data represents a different population from the old data: while the two have different bins, the new data suggests higher success for most data. This may be from a different method of collecting the data, or golfers improving in the intervening years.

Fitting the model on the new data

Since we think these may be two different populations, the easiest solution would be to refit our model. This goes worse than earlier: there are divergences, and it takes much longer to run. This may indicate a problem with the model: Andrew Gelman calls this the "folk theorem of statistical computing".

In [19]:
with angle_model(new_golf_data):
    new_angle_trace = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [variance_of_shot]
100.00% [8000/8000 00:15<00:00 Sampling 4 chains, 186 divergences]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 16 seconds.
There were 11 divergences after tuning. Increase `target_accept` or reparameterize.
There were 143 divergences after tuning. Increase `target_accept` or reparameterize.
There were 13 divergences after tuning. Increase `target_accept` or reparameterize.
There were 19 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.
In [20]:
ax = plot_golf_data(new_golf_data)
plot_golf_data(golf_data, ax=ax)
t = np.linspace(0, new_golf_data.distance.max(), 200)

ax.plot(t, forward_angle_model(angle_trace['variance_of_shot'].mean(), t), label='Trained on original data')
ax.plot(t, forward_angle_model(new_angle_trace['variance_of_shot'].mean(), t), label='Trained on new data')
ax.set_title('Retraining the model on new data')
ax.legend();

A model incorporating distance to hole

We might assume that, in addition to putting in the right direction, a golfer may need to hit the ball the right distance. Specifically, we assume:

  1. If a put goes short or more than 3 feet past the hole, it will not go in.
  2. Golfers aim for 1 foot past the hole
  3. The distance the ball goes, $u$, is distributed according to $$ u \sim \mathcal{N}\left(1 + \text{distance}, \sigma_{\text{distance}} (1 + \text{distance})\right), $$ where we will learn $\sigma_{\text{distance}}$.

Again, this is a geometry and algebra problem to work the probability that the ball goes in from any given distance: $$ P(\text{good distance}) = P(\text{distance} < u < \text{distance} + 3) $$

it uses Phi, the cumulative normal density function we implemented earlier.

In [21]:
OVERSHOT = 1.0
DISTANCE_TOLERANCE = 3.0

def distance_angle_model(golf_data):
    distances = golf_data.distance.values
    with pm.Model() as distance_angle_model:
        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(tt.arcsin((CUP_RADIUS - BALL_RADIUS) / distances) / variance_of_shot) - 1)
        p_good_distance = pm.Deterministic('p_good_distance', Phi((DISTANCE_TOLERANCE - OVERSHOT) / ((distances + OVERSHOT) * variance_of_distance)) - 
                                                              Phi(-OVERSHOT / ((distances + OVERSHOT) * variance_of_distance)))

        success = pm.Binomial('success', n=golf_data.tries, p=p_good_angle * p_good_distance, observed=golf_data.successes)
    return distance_angle_model

pm.model_to_graphviz(distance_angle_model(new_golf_data))
Out[21]:
%3 cluster31 31 variance_of_distance variance_of_distance ~ HalfNormal p_good_distance p_good_distance ~ Deterministic variance_of_distance->p_good_distance variance_of_shot variance_of_shot ~ HalfNormal p_good_angle p_good_angle ~ Deterministic variance_of_shot->p_good_angle success success ~ Binomial p_good_distance->success p_good_angle->success

This model still has only 2 dimensions to fit. We might think about checking on OVERSHOT and DISTANCE_TOLERANCE. Checking the first might involve a call to a local golf course, and the second might require a trip to a green and some time experimenting. We might also think about adding some explicit correlations: it is plausible that less control over angle would correspond to less control over distance, or that longer putts lead to more variance in the angle.

Fitting the distance angle model

In [22]:
with distance_angle_model(new_golf_data):
    distance_angle_trace = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [variance_of_distance, variance_of_shot]
100.00% [8000/8000 00:05<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 6 seconds.
The acceptance probability does not match the target. It is 0.8798848192814278, but should be close to 0.8. Try to increase the number of tuning steps.
The acceptance probability does not match the target. It is 0.8892228299255789, but should be close to 0.8. Try to increase the number of tuning steps.
The number of effective samples is smaller than 25% for some parameters.
In [23]:
def forward_distance_angle_model(variance_of_shot, variance_of_distance, t):
    rv = 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)
t = np.linspace(0, new_golf_data.distance.max(), 200)


ax.plot(t, forward_angle_model(new_angle_trace['variance_of_shot'].mean(), t), label='Just angle')
ax.plot(t, forward_distance_angle_model(distance_angle_trace['variance_of_shot'].mean(), distance_angle_trace['variance_of_distance'].mean(), t), label="Distance and angle")

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

This new model looks better, and fit much more quickly with fewer sampling problems compared to the old model.There is some mismatch between 10 and 40 feet, but it seems generally good. We can come to this same conclusion by taking posterior predictive samples, and looking at the residuals. Here, we see that the fit is being driven by the first 4 bins, which contain ~40% of the data.

In [24]:
with distance_angle_model(new_golf_data):
    ppc = pm.sample_posterior_predictive(distance_angle_trace)
    
residuals = 100 * (new_golf_data.successes - ppc['success'].mean(axis=0)) / new_golf_data.tries

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:07<00:00]

A new model

It is reasonable to stop at this point, but if we want to improve the fit everywhere, we may want to choose a different likelihood from the Binomial, which cares deeply about those points with many observations. One thing we could do is add some independent extra error to each data point. We could do this in a few ways:

  1. The Binomial distribution in usually parametrized by $n$, the number of observations, and $p$, the probability of an individual success. We could instead parametrize it by mean ($np$) and variance ($np(1-p)$), and add error independent of $n$ to the likelihood.
  2. Use a BetaBinomial distribution, though the error there would still be (roughly) proportional to the number observations
  3. Approximate the Binomial with a Normal distribution of the probability of success. This is actually equivalent to the first approach, but does not require a custom distribution. Note that we will use $p$ as the mean, and $p(1-p) / n$ as the variance. Once we add some dispersion $\epsilon$, the variance becomes $p(1-p)/n + \epsilon$.

We follow approach 3, as in the Stan case study, and leave 1 as an exercise.

In [25]:
def disp_distance_angle_model(golf_data):
    distances = golf_data.distance.values
    with pm.Model() as distance_angle_model:
        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(tt.arcsin((CUP_RADIUS - BALL_RADIUS) / distances) / variance_of_shot) - 1)
        p_good_distance = pm.Deterministic('p_good_distance', Phi((DISTANCE_TOLERANCE - OVERSHOT) / ((distances + OVERSHOT) * variance_of_distance)) - 
                                                              Phi(-OVERSHOT / ((distances + OVERSHOT) * variance_of_distance)))
        
        p = p_good_angle * p_good_distance
        p_success = pm.Normal('p_success', mu=p, sd=tt.sqrt(((p * (1 - p)) / golf_data.tries) + dispersion ** 2), observed=golf_data.successes / golf_data.tries)
    return distance_angle_model

pm.model_to_graphviz(disp_distance_angle_model(new_golf_data))
Out[25]:
%3 cluster31 31 variance_of_distance variance_of_distance ~ HalfNormal p_good_distance p_good_distance ~ Deterministic variance_of_distance->p_good_distance dispersion dispersion ~ HalfNormal p_success p_success ~ Normal dispersion->p_success variance_of_shot variance_of_shot ~ HalfNormal p_good_angle p_good_angle ~ Deterministic variance_of_shot->p_good_angle p_good_distance->p_success p_good_angle->p_success
In [26]:
with disp_distance_angle_model(new_golf_data):
    disp_distance_angle_trace = pm.sample(1000, tune=1000)
    disp_ppc = pm.sample_posterior_predictive(disp_distance_angle_trace)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [dispersion, variance_of_distance, variance_of_shot]
100.00% [8000/8000 00:05<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 6 seconds.
100.00% [4000/4000 00:14<00:00]
In [27]:
ax = plot_golf_data(new_golf_data, ax=None)
t = np.linspace(0, new_golf_data.distance.max(), 200)


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