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.