Bayesian Estimators for the Beta-Binomial Model of Batting Ability

In my last post introducing Bayesian stats through the simplest non-trivial distribution, the binomial, I introduced moment-matching “empirical Bayes” point estimates of the beta priors. In this post, I’ll introduce the so-called “Bayesian estimator” point estimate for the beta priors. In addition, I’ll show why maximum a posteriori (MAP) estimates of batting average differ from the Bayesian estimates for batting average.

Bayesian Estimators

Suppose we provide an estimate $\hat{\theta}$ for a parameter that has true value $\theta$ . We define a loss function $\mathcal{L}(\hat{\theta},\theta)$ , lower values of which indicate that $\hat{\theta}$ is a better estimate of $\theta$ .

The Bayesian estimator minimizes the expected loss given a Bayesian posterior distribution $p(\theta|x)$ over parameter $\theta$ :

$\theta^{\mathcal{L}} = \arg\min_{\theta} \mathbb{E}_{p(\theta|x)}[\mathcal{L}(\hat{\theta},\theta)]$

$= \arg\min_{\theta} \int_{\Theta} \mathcal{L}(\hat{\theta},\theta) \, p(\theta|x) \, d\theta$

With Squared Loss, Bayesian Estimates are Posterior Means

The most common loss function, due to its convenience, is squared loss:

$\mathcal{L}(\hat{\theta},\theta) = L_2(\theta,\hat{\theta}) = (\theta - \hat{\theta})^2$

As stated, this assumes a single real-valued parameter; the natural extension is to square distance for parameter vectors.

With squared loss, the expected loss is minimized by the mean of the posterior. In other words, the Bayes estimator for a parameter under squared loss is the posterior mean of that parameter:

$\theta^{L_2} = \int_{\Theta} \theta \, p(\theta|x) \, d\theta$

Beta-Binomial Batting Model

Our model for batting so far is very simple, with player $j$ ‘s ability being drawn from a beta prior with fixed hyperparameters $\alpha$ (prior hits plus 1) and $\beta$ (prior outs plus 1):

$\theta_j \sim \mbox{\sf Beta}(\alpha,\beta)$

The number of hits $n_j$ for player $j$ in $N_j$ at bats is drawn from a binomial sampling distribution:

$n_j \sim \mbox{\sf Binom}(\theta_j)$

The observed batting average is just $n_j/N_j$ .

Recall that the binomial is defined by:

$\mbox{\sf Binom}(n|N,\theta) \propto \theta^n (1-\theta)^{N-n}$

and the beta prior by:

$\mbox{\sf Beta}(\theta|\alpha,\beta) \propto \theta^{\alpha} (1 - \theta)^{\beta}$

Given the convenient exponential form of the prior and likelihood function, the posterior is easily computed as:

$p(\theta_j|n_j,N_j,\alpha,\beta) \propto p(n_j|N_j,\theta_j) \, p(\theta_j|\alpha,\beta)$

$= \mbox{\sf Binom}(n_j|N_j,\theta_j) \, \mbox{\sf Beta}(\theta_j|\alpha,\beta)$

$\propto (\theta_j^{n_j}(1-\theta_j)^{N_j - n_j}) \, (\theta_j^{\alpha} (1 - \theta_j)^{\beta}$

$= \theta_j^{n_j + \alpha} \, (1 - \theta_j)^{N_j - n_j + \beta}$

Given this final form, which is proportional to a beta, we have

$\theta_j \sim \mbox{\sf Beta}(n_j+\alpha, N_j-n_j+\beta)$ .

Thus the beta prior is conjugate to the binomial distribution, because when the prior over abilities is a beta density and the sampling distribution for the number of hits a binomial distribution, the posterior over abilities is also a beta density.

Bayesian Estimator for Batting Ability

The posterior mean for batting ability $\theta_j$ for player $j$ is thus the mean of the posterior beta distribution. The mean of a beta distribution $\mbox{\sf Beta}(\alpha,\beta)$ is $\frac{\alpha}{\alpha+\beta}$ . So for the posterior $\mbox{\sf Beta}(n_j+\alpha, N_j-n_j+\beta)$ , the Bayes estimators for batting ability are just:

$\theta_j^{L_2} = \frac{n_j + \alpha}{N_j + \alpha + \beta}$

Contrast with the MAP Estimate

Note that this is not the same estimate as the maximum a posteriori (MAP) estimate, because the mode (maximum point) and mean of a beta distribution are not the same. The MAP estimate is based on the mode,

$\theta_j^* = \frac{n_j + \alpha-1}{N_j + \alpha + \beta-2}$

Diffuse Priors for the Beta Priors

In order to derive a Bayes estimate of the prior parameters $\alpha$ and $\beta$ , we need to extend our model with a prior for the priors (as we said last time, it’s elephants all the way down). As is common in these priors-for-priors situations, we’ll put a very diffuse prior on the priors themselves.

For convenience, we’ll reparameterize the beta distribution using its mean $\alpha/(\alpha+\beta)$ and scale $\alpha + \beta$ . Given the mean and scale, we can derive $\alpha$ and $\beta$ by:

$\alpha = (\alpha + \beta) (\frac{\alpha}{\alpha+\beta})$

$\beta = (\alpha + \beta) (1 - \frac{\alpha}{\alpha+\beta})$ .

We then put priors directly on the mean and scale. Because the mean has a range of 0 to 1, we can use a uniform prior:

$\frac{\alpha}{\alpha+\beta} \sim \mbox{\sf Uniform}(0,1) = \mbox{\sf Beta}(1,1)$

Because the scale ranges from 0 to infinity, it’s impossible to have a (proper) uniform prior, because the integral would diverge. Instead, we use a Pareto prior,

$(\alpha + \beta) \sim \mbox{\sf Pareto}(1.5,1)$

The Pareto distribution $\mbox{\sf Pareto}(\gamma,k)$ is defined on the interval $[k,\infty)$ by

$\mbox{Pareto}(y|\gamma,k) \sim y^{-(\gamma+1)}$

In other words, the Pareto distribution is a power-law distribution with exponent $\gamma+1$ . This is a fairly diffuse prior, slightly favoring smaller values of the scale.

BUGS Code

Here’s the full model coded in BUGS:

model {
  for (j in 1:J) {
    ability[j] ~ dbeta(prior.hits,prior.outs)
    hits[j] ~ dbin(ability[j],atbats[j])
  }
  prior.hits <- prior.avg * prior.atbats
  prior.outs <- (1 - prior.avg) * prior.atbats
  prior.avg ~ dbeta(1,1)
  prior.atbats ~ dpar(1.5,1)
}

This model assumes J batters, with abilities $\theta_j$ coded as ability[j], hits $n_j$ as hits[j], at bats $N_j$ as atbats[j]. The prior parameter $\alpha$ is coded as prior.hits, and $\beta$ as prior.outs. The reparameterized prior mean $\alpha/(\alpha+\beta)$ is coded as prior.avg, and the reparameterized prior scale $\alpha+\beta$ as prior.atbats.

Gibbs Sampling

With the model coded in BUGS, we can automatically run Gibbs sampling to collect a number of samples from the posterior distribution $p(\alpha,\beta,\theta|n,N)$ (we suppressed the conditioning on the hyperparameters).

We can then approximate the posterior marginal averages (to any number of decimal places) by taking the average of the samples for the respective variables.

It’s very easy to call BUGS from R using the >R2WinBUGS package. So we read the data into R, use it to call BUGS, then we collect the results back in R, where we can calculate their means.

data.tsv = scan("2006ALnonPitch.tsv")
N = length(data.tsv)/2
y = matrix(data.tsv, nrow=N, ncol=2, byrow=TRUE)
hits = y[,1]
atbats = y[,2]

library("R2WinBUGS")
data = list("atbats","hits","N")
parameters = c("ability","prior.avg","prior.atbats")

inits = function() {
  list(ability=rep(N,0.30),
       prior.hits=c(1),prior.outs=c(2),
       prior.avg=c(1.0/3.0), prior.atbats=c(3))
}
anno = bugs(data, inits, parameters,
    "c:/carp/devmitz/carp/baseball/src/baseballone/binomial.bug",
    n.chains=3, n.iter=1000,
    debug=TRUE,
    clearWD=TRUE,
    bugs.directory="c:\WinBUGS\WinBUGS14")

print(anno)
plot(anno)
attach.bugs(anno)

This presupposes that you’ve installed BUGS at the location specified by bugs.directory, have provided the BUGS program at the specified full path (I can’t get relative paths to work here), and have provided the data in the file 2006ALnonPitch.tsv (see the last section for a dump of the data in this file). The format of the tsv (tab-separated values) file is one player per line, with the number of hits, a tab, the number of at bats, and a newline.

The parameters specify three simultaneous Markov chains for 1000 iterations each. By default, BUGS will discard the first half of the iterations to remove noise from before the chains converge.

We initialize the ability values to $\theta_j = 0.300$ , $\alpha=1$ and $\beta=2$ . Better initialization helps models converge to the posterior samples more quickly, but aren’t even necessary for easy-to-fit models like these that run relatively quickly.

The data and parameters are specified. Data must be defined in R before calling BUGS, and parameters are returned as a two-dimensional array of samples (indexed by chain and sample number).

BUGS is not super-fast, taking about 15 minutes to take 3000 samples in 3 chains on my old notebook (it’s about twice as fast on my newish workstation).

Bayesian Estimates

While I don’t have time to explain Gibbs sampling in detail in this context, let me just say that the three chains mixed very well with $\hat{R}$ values all having values very near 1.0.

The mean prior batting average is 0.271, and the mean prior at bats was 403. Note that this average is lower and the number of at bats higher than with the moment-matching estimates. I calculated these values in R after closing the BUGS window by:

> mean(prior.avg)
[1] 0.2709107
> mean(prior.atbats)
[1] 403.1252

This corresponds to beta density parameters $\alpha=109.2$ and $\beta=293.9$ .

I’ll discuss more about the joint marginal average/at-bats posterior and the marginal ability posteriors in the next and (hopefully) final post in the series.

Who Do We Think is Better?

Who do we think is a better batter, someone we’ve seen go 40/100 (a 0.400 average) or someone we’ve seen go 175/500 (a 0.350 average)?

For now, suppose we fix an empirical Bayes prior to the Bayesian estimate $\mbox{\sf Beta}(109.2,293.9)$ , so that $\alpha=109.2$ and $\beta=293.9$ The Bayesian estimator (posterior mean) for a player’s ability $\theta$ when we observe the player getting get $n$ hits in $N$ at bats, is then the expected value of the posterior, which is:

$\theta^{L_2} = \frac{n+\alpha-1}{N + \alpha + \beta - 2}$

If we plug in $n=40, N=100$ we get an estimated ability of 0.295. If we plug in $n=175, N=500$ , we get an estimated ability of 0.314! If we believe this model, we can’t just compare batting averages — the number of at bats is a critical element of our posterior estimate for low numbers of at bats. This is the effect of the prior, not the Bayesian inference per se. Of course, the scale of the prior being so large is the result of using the Bayesian estimator for the prior rather than the moment-matching estimate, so the effect is greater.

In the Limit

In the limit, of course, the Bayesian estimate approaches the maximum likelihood estimate:

$\lim_{N \rightarrow \infty} \frac{n + \alpha + 1}{N + \alpha + \beta + 1} = \frac{n}{N}$

Data from Baseball1

I gathered the data for the 2006 American League position players from Baseball1 (click on the “download database” menu item on the nav bar). It requires merging the player table (to get positions and exclude piters) and the batting table. For completeness, here’s the munged data:

This entry was posted on September 23, 2009 at 11:08 am and is filed under Carp's Blog, Statistics. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

LingPipe Blog