Choose a prior’s parameters

How to parameterise a prior when all you know is the mean and one other percentile.

Duncan Garmonsway
December 14, 2018

Updated 2019-03-05: Python implementation and R improvements.

Given a distribution, how do you find the parameters? When choosing a beta prior for a Bayesian A/B test, for example.

I came up with these parameters for two reasons:

• The mean is $$\frac{\alpha}{\alpha + \beta} = \frac{81}{81 + 219} = .270$$

• As you can see in the plot, this distribution lies almost entirely within (.2,.35)- the reasonable range for a batting average. – David Robinson

What David Robinson is doing is best-practice for choosing a prior – ask an expert. But he makes a huge leap from what the expert says to what the parameters are. Normally you’d use dbeta(shape1, shape2) to get the distribution for given parameters, but here the problem is reversed: find the parameters for a distribution with a mean of .270 and a cumulative density of, say, .99 between .2 and .35.

How did David Robinson come up with the parameters 81 and 219? There aren’t tools for that. Perhaps he’s brainy enough to do the maths in his head and think, “Hmm, beta with mean = .270 and most of it between .2 and .35, in that case 81 and 219 should be about right.” Or perhaps he just tried some numbers and adjusted them until they worked. There are online tools for that, or you can use R. For example, here’s a beta distribution with an initial guess of the parameters shape1 = 2 and shape2 = 4).

curve(dbeta(x, 2, 4), from = 0, to = 1, type = "l") Not quite right – the parameters need adjusting. And it turns out to be fairly straightforward to find parameters that get exactly the right mean by working with the formula for the mean $$\frac{\alpha}{\alpha + \beta}$$. The mean can be interpreted as a number of successes divided by a number of attempts (successes + failures), so a mean of .270 successes every 1 attempt is like 270 successes every 1000 attempts. R’s parameterisation of the beta distribution is shape1 successes and shape2 failures, so a simple subtraction converts 270 successes and 100 attempts into (1000 - 270) failures.

curve(dbeta(x, 270, (1000 - 270)), from = 0, to = 1, type = "l") Now the distribution has correct mean, but the peak is too narrow – most of it lies above 0.25, whereas David Robinson only aimed for it to be above .2 and below .35. Unfortunately the cumulative density function is to complicated to hack, but using a handy property – smaller parameters make wider peaks – let’s divide the parameters by a factor of, say, 3.

curve(dbeta(x, 270 / 3, (1000 - 270) / 3), from = 0, to = 1, type = "l") 270 / 3 is 90 for $$\alpha$$, and (1000 - 270) / 3 is about 243 for $$\beta$$, so this is in the ballpark of David Robinson’s parameters 81 and 219, and we could stop here.

No don’t stop! Don’t you see? What just happened was tedious trial and error; manual optimisation of the parameters to get the desired shape! The computer used you as a labour-saving device. Give a taste of its own medicine! Make the computer do the trial and error with optim()!

Here’s a function that plugs various values of shape1 and shape2 into qbeta() until it finds a curve with the properties we want: a mean of .270 and, say, 1% of the distribution to be below .2.

Effectively this re-parameterises the beta distribution to use “mean” and “amount below x” instead of shape1 “successes” and shape2 “failures.” Call it the Gimme Beta for ‘Gimme the beta I want’. But because I am lazy have a computer, instead of sweating over algebra to translate “mean” and “amount below x” into shape1 and shape2, I use optim() to find values of shape1 and shape2 that create a distribution that also satisifies my “mean” and “amount below x”. And then I can use shape1 and shape2 to do the A/B test, or whatever the prior is for.

$\mathrm{Gimme~Beta}(\mathrm{mean}, \mathrm{amount~below~}x) \sim \mathrm{Beta}(\mathrm{shape1}, \mathrm{shape2})$

Updated 2019-03-05: Use absolute differences rather than squared differences, because squared differences shrink to small for floating-point arithmetic to be accurate.

gimme_beta <- function(mean, cumulative_density, below, par_start = c(2, 2),
plot_from_to = NULL) {
error <- function(x) {
alpha1 <- x
beta1  <- x
fitted_mean <- alpha1 / (alpha1 + beta1)
fitted_cumulative_density <- pbeta(below, alpha1, beta1)
diff_mean <- abs(fitted_mean - mean)
diff_cumulative_density <- abs(fitted_cumulative_density - cumulative_density)
diffs <- sum(diff_mean,
diff_cumulative_density)
diffs
}
estimate <- optim(par_start, error)
alpha1 = estimate$par beta1 = estimate$par
if (is.null(plot_from_to)) {
plot_from_to <- c(qbeta(0.00005, alpha1, beta1),
qbeta(0.99995, alpha1, beta1))
}
y <- curve(dbeta(x, alpha1, beta1),
from = plot_from_to,
to = plot_from_to,
type = "l")
c(shape1 = alpha1, shape2 = beta1)
}

How well does gimme_beta() do on David Robinson’s task?

gimme_beta(mean = .270, cumulative_density = .01, below = .2) shape1    shape2
52.75807 142.64219

That isn’t far off the 81 and 219 that David Robinson came up with, but then we don’t know exactly what David Robinson meant by “lies almost entirely within”. What if we work out the $$\mu$$ and $$\theta$$ parameters of David Robinson’s actual distribution, plug them back into gimme_beta() and check again.

81 / (81 + 219) # mu, the mean

 0.27

pbeta(.2, 81, 219) # theta, the percentage of the distribution that is below 0.2

 0.001897932

So plugging the (kinda) mean and the percentage-that-is-below-0.2 into gimme_beta(), does it give us David Robinson’s parameters back?

gimme_beta(mean = .270, cumulative_density = .001897932, below = .2) shape1 shape2
81    219

Close enough for government work. And more importantly, gimme_beta() will gimme the parameters I need for my A/B test, which happen to be a mean of .23 with 5% of the distribution below 0.2.

gimme_beta(.23, .05, .2) shape1   shape2
117.3925 393.0097

Updated 2019-03-05: Python implementation.

Thanks to Sarah Dean for pair-coding this translation.

from scipy.stats import beta
from scipy.optimize import minimize

def f(x, desired_mean, cumulative_density, below_point):
fitted = beta(x, x)
fitted_mean = fitted.mean()
fitted_cdf = fitted.cdf(below_point)
diff_mean = abs(fitted_mean - desired_mean)
diff_cdf = abs(fitted_cdf - cumulative_density)
diff = diff_mean + diff_cdf
return diff

fitted = minimize(f, [2.0, 2.0], args = (0.270, 0.01, 0.2), method = 'Nelder-Mead')
fitted_dist = beta(fitted.x, fitted.x)

print(fitted.x)

[ 52.75811484 142.64234025]

print(fitted_dist.mean())

0.26999995889124706

print(fitted_dist.cdf(0.2))

0.009999998660408616

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/nacnudus/duncangarmonsway, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

Garmonsway (2018, Dec. 14). Duncan Garmonsway: Choose a prior's parameters. Retrieved from https://nacnudus.github.io/duncangarmonsway/posts/2018-12-14-choose-a-priors-parameters/

BibTeX citation

@misc{garmonsway2018choose,
author = {Garmonsway, Duncan},
title = {Duncan Garmonsway: Choose a prior's parameters},
url = {https://nacnudus.github.io/duncangarmonsway/posts/2018-12-14-choose-a-priors-parameters/},
year = {2018}
}