# Is Bayesian model selection with empirical parameter priors sound?

### 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.

• A coin with $\theta=0.52$ counts as loaded, doesn't it? "Ignoring coins with $0.4<\theta<0.6$" will likely ignore a lot of loaded coins (at least if there are enough overall), and will probably push you to overestimate the amount of bias of the loaded coin, and to miss mildly loaded coins. In any case, "Ignoring coins with $0.4<\theta<0.6$" assumes that you know $\theta$, which you don't. Chances are you mean some estimator $\hat\theta$. Apr 15, 2023 at 22:35
• @ChristianHennig: Ignoring coins with 0.4<𝜃<0.6" will likely … – Mind that this is just to estimate the form of the prior, which I can (and must) assume to be reasonably continued for $0.4<θ<0.6$ (as opposed to dropping to zero in this region). I could also fit a mixture distribution such as described in J. Delaney’s answer, but that would require a good a priori estimate of $p_\text{fair}$. As for your second point, please see my edit. Apr 16, 2023 at 0:12
• Is the idea that there is a fixed probability for a loaded coin and the machine uses this probability whenever it produces a loaded coin, or does the machine produce loaded coins with (randomly?) varying probabilities? In other words, does the loaded prior for $\theta$ model an empirical distribution of $\theta$-values produced by the machine, or is it rather an epistemic prior for a single value of $\theta$ that is actually realised? Apr 16, 2023 at 10:55
• @ChristianHennig: The latter. I state this more explicitly now. Apr 16, 2023 at 11:45

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.

• Thank you for your answer. I am not entirely sure whether I get the main direction of it: Are you saying that 1) what you describe is essentially equivalent to my approach and I just have miscategorised it or 2) what you describe is a different approach? In the latter case, it seems to me that the two approaches are equivalent at first glance (except maybe for different ways to point estimate) and the model-selection approach has the advantage that it allows a certain dissection of the problem that is convenient for more complex cases. If not, what would be bad about my approach? Apr 15, 2023 at 21:18
• It's the former (1). Your approach is similar to an empirical Bayes method, except that the method of estimation can be made more rigorous (no need to arbitrarily ignore some of the coins). Apr 16, 2023 at 8:53
• Thanks. I am still somewhat skeptical of the approach with the mixture distribution being generally better for the following reason: If the distribution for loaded coins has a prominent peak at ½ (thus violating my constraint of being nicely smooth) this is difficult to distinguish from a higher prior for fair coins. For this to work, I would need a very good prior for fair coins or a lot of very accurate data (neither of which I have in my actual application). Apr 16, 2023 at 18:03
• If you are worried about the accuracy of the point estimates due to insufficient data, then you should take the fully Bayesian approach (the method of ignoring the central estimates will suffer from the same issue). Notice that you can estimate the uncertainty of the point estimates by analyzing the likelihood function. If you require symmetry by imposing $\alpha=\beta$, then you are left with a only 2-parameter model, so you can even visualize the contours of $\mathcal L(p_f,\alpha)$. It is also easy enough to run simulations and compare the estimates to the true values. Apr 16, 2023 at 18:31