Is Bayesian model selection with empirical parameter priors sound?

Question

Overview

I want to perform a Bayesian model selection on many datasets and use these datasets to determine the required parameter priors.

Example Scenario: Coins

Suppose I have a collection of thousand coins produced by a machine that randomly produces fair and loaded coins. The loaded coins are not identical, but their heads ratio $θ$ follows an unknown distribution $p(θ|\mathcal{M}_\text{loaded})$ obeying some constraints (see below). For each coin, I want to decide whether it’s fair using Bayesian model selection with two models $\mathcal{M}_\text{loaded}$ and $\mathcal{M}_\text{fair}$.

I know:

• For each coin: the number of heads from hundred tosses (and thus an estimator $\hat{θ}$ for the heads ratio $θ$).
• Model priors $p_\text{fair}$ and $p_\text{loaded}$ with $0.1≤p_\text{fair}≤0.9$.
• The probability density $p(θ|\mathcal{M}_\text{loaded})$ of the heads ratio of the loaded coins obeys the following constraints:
  • symmetric around ½,
  • smooth,
  • not very far from a uniform distribution, say, $0.1 < p(θ|\mathcal{M}_\text{loaded}) < 10$ everywhere.

With all this given, the main information I am lacking for this is a prior for $p(θ|\mathcal{M}_\text{loaded})$. I estimate this by finding a suitable distribution and fitting it to my data for all coins, ignoring coins with $0.4<\hat{θ}<0.6$, since those have a decent chance to be fair. The rest of the Bayesian model selection is straightforward.

Questions

• Is this procedure sound? I acknowledge that I use the same data twice. However, the data for a given coin has barely any impact on the parameter priors relevant to its model selection. (I could also exclude the data for the given coin when determining the priors for its analysis, doing a thousand fits instead of just one.)

• If yes, is there a name or reference for this approach?

• If no, is there a better way to determine parameter priors for $\mathcal{M}_\text{loaded}$? I am particularly interested in ways that can be extended to a more complex model space as well as higher-dimensional and unbounded parameter spaces.

Answer 1

The scenario in your example is actually better described not as "model selection" problem (where you have to decide between two models to describe the entire data) but rather as an Empirical Bayes method applied to a hierarchical model.

Specifically you assume that $\theta$ has a mixture distribution

$$ \theta \sim p_f\delta(1/2) + (1-p_f)\text{Beta}(\theta|\alpha,\beta)$$

if for example you choose a Beta distribution to describe $p(\theta|\mathcal M_{loaded})$. Then you can use this to estimate the model parameters $p_f,\alpha,\beta$ by calculating the marginal likelihood:

$$P(\text{data}| p_f,\alpha,\beta) = \prod_i \int d\theta_i P(\text{data}|\theta_i)\times P(\theta_i | p_f,\alpha,\beta). $$

The "empirical" aspect in "empirical Bayes" refers to using point estimates of those parameters (for example by maximum likelihood) as a prior for a particular $\theta_i$ (the "low" level in the hierarchy), for example if you have count data $k_i \sim \text{Binomial}(n_i,\theta_i)$ then you would calculate the posterior probability of $\theta_i$ as

$$ P(\theta_i|k_i,n_i) \propto \theta_i^{k_i} (1-\theta_i)^{n_i-k_i} P(\theta_i | \hat p_f,\hat \alpha,\hat \beta)$$

where a "hat" over a parameter denotes its point estimate${^1}$.

It is also possible to treat this in a "fully Bayesian" way by assigning a prior $\pi(p_f,\alpha,\beta)$ to the high level parameters and marginalizing, using the full posterior distribution:

$$ P(\theta_i|\text{data}) \propto \theta_i^{k_i} (1-\theta_i)^{n_i-k_i} \int d\alpha\int d\beta \int dp_f P(\theta_i|p_f,\alpha,\beta)\times P(p_f,\alpha,\beta | \text{data}_{-i}) $$

Where

$$P(p_f,\alpha,\beta | \text{data}_{-i}) \propto \int d\theta_1 ... \int d\theta_n \prod_{j\ne i} \theta_j^{k_j} (1-\theta_j)^{n_j-k_j} P(\theta_j|p_f,\alpha,\beta)\pi(p_f,\alpha,\beta)$$

Is the posterior distribution of the hyperparameters. To be completely rigorous indeed requires excluding the coin of interest from the data, however for large dataset this might be a negligible effect. Calculating the integrals in hierarchical models can usually be done only numerically (for example the Beta distribution does not have a simple conjugate prior). However when the dataset is large enough we can expect the posterior probability to become concentrated around the point estimates, such that the full calculation reduces to the empirical one. This is the justification for using empirical Bayes methods as an approximation.

---

${^1}$ Using the Beta distribution is convenient since it is the conjugate prior of the Binomial distribution, so the marginal likelihood can be calculated analytically:

$$ \int d\theta \theta^k (1-\theta)^{n-k}\times (p_f\delta(1/2) + (1-p_f)\text{Beta}(\theta|\alpha,\beta)) $$

$$=p_f\frac{1}{2^n} + (1-p_f)\frac{B(\alpha+k,\beta+n-k)}{B(\alpha,\beta)}$$

where $B(\cdot,\cdot)$ is the Beta function. The maximum likelihood estimates are

$$\hat p_f,\hat \alpha,\hat \beta = \underset{p_f,\alpha,\beta}{\text{argmax}} \sum_i \log \left( p_f\frac{1}{2^{n_i}} + (1-p_f)\frac{B(\alpha+k_i,\beta+n_i-k_i)}{B(\alpha,\beta)} \right)$$

and the posterior distribution of $\theta_i$ is

$$P(\theta_i|k_i,n_i) = \tilde p_f \delta(1/2) + (1-\tilde p_f)\text{Beta}(\theta_i|\hat \alpha+k_i, \hat \beta+n_i-k_i)$$

where

$$\tilde p_f = \frac{\hat p_f\frac{1}{2^{n_i}}}{ \hat p_f\frac{1}{2^{n_i}} + (1-\hat p_f)\frac{B(\hat \alpha+k_i,\hat \beta+n_i-k_i)}{B(\hat \alpha,\hat \beta)}} $$

is the posterior probability of the coin being fair.