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) + p_l\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,p_l,\alpha,\beta$ by calculating the marginal likelihood:
$$P(\text{data}| p_f,p_l,\alpha,\beta) = \prod_i \int d\theta_i P(\text{data}|\theta_i)\times P(\theta_i | p_f,p_l,\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) \propto \theta_i^{k_i} (1-\theta_i)^{n_i-k_i} P(\theta_i | \hat p_f,\hat p_l,\hat \alpha,\hat \beta)$$
where a "hat" over a parameter denotes its point estimate.
It is also possible to treat this in a "fully Bayesian" way by using an additional prior for the high level parameters $\pi(p_f,p_l,\alpha,\beta)$ and marginalizing:
$$ 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\int dp_l P(\theta_i|p_f,p_l,\alpha,\beta)\times P(p_f,p_l,\alpha,\beta | \text{data}_{-i}) $$
Where
$$P(p_f,p_l,\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,p_l,\alpha,\beta)\pi(p_f,p_l,\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.