Additional-Review-(marginalization-over-joint-population+EoS-uncertainty).md

@@ -32,14 +32,42 @@ The new version of the code relies on the user to specify the draw probability,
 
 **bad population reweighing**
 
-The previous version of the code implemented population priors by reweighing single-event PE samples. This was done within [this code block](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/c8b0d297d11b0a350c23ee44832d24b899729d8f/bin/mmax-model-selection#L227-270). Importantly, we see that the total weight for each PE sample was computed as the sum of population priors (implicit from [this block](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/c8b0d297d11b0a350c23ee44832d24b899729d8f/bin/mmax-model-selection#L263-270)). That sum was then normalized by the default PE prior [here](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/c8b0d297d11b0a350c23ee44832d24b899729d8f/bin/mmax-model-selection#L279).
+The previous version of the code implemented population priors by reweighing single-event PE samples. This was done within [this code block](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/c8b0d297d11b0a350c23ee44832d24b899729d8f/bin/mmax-model-selection#L227-270). Importantly, we see that the total weight for each PE sample was computed as the sum of population priors (implicit from [this block](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/c8b0d297d11b0a350c23ee44832d24b899729d8f/bin/mmax-model-selection#L263-270)). That sum was then normalized by the default PE prior [here](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/c8b0d297d11b0a350c23ee44832d24b899729d8f/bin/mmax-model-selection#L279). This means the overall weight assigned to each PE sample was
 
-However, as is shown in the [technical note](#technical
+```math
+w_i = \frac{\sum_p p(\theta_i)|\Lambda_p)}{p(\theta_i|\mathrm{default\ PE})}
+```
+
+This implies that the total sum would be
+
+```math
+S_\mathrm{old} = \sum_i w_i \Theta_i = \sum_i \left(\frac{\sum_p p(\theta_i)|\Lambda_p)}{p(\theta_i|\mathrm{default\ PE})}\right) \Theta_i = \sum_p \sum_i \frac{p(\theta_i|\Lambda_p)}{p(\theta_i|\mathrm{default\ PE})} \Theta_i
+```
+
+It is straightforward to show that
+
+```math
+\sum_i \frac{p(\theta_i|\Lambda)}{p(\theta_i|\mathrm{default\ PE})} \Theta_i \approx \int d\theta p(\theta|\mathrm{data},\Lambda) \Theta(\theta) \frac{p(\mathrm{data}|\Lambda)}{p(\mathrm{data}|\mathrm{default\ PE})}
+```
+
+which implies the total sum would approximate
+
+```math
+S_\mathrm{old} \approx \int d\Lambda p(\Lambda) \frac{p(\mathrm{data}|\Lambda)}{p(\mathrm{data}|\mathrm{default\ PE})} \int d\theta p(\theta|\mathrm{data},\Lambda) \Theta(\theta)
+```
+
+and there is an extra factor of `p(data|Lambda)` included within the integrand.
+As is shown in the [technical note](#technical-note), we instead wish to compute
+
+```math
+S = \int d\Lambda p(\Lambda) \int d\theta p(\theta|\mathrm{data},\Lambda) \Theta(\theta)
+```
 
+and to do so we must normalize the sum over `\theta_i` to remove the factor of `p(data|Lambda)`.
+The correct sum is implemented within the updated code [here](https://git.ligo.org/reed.essick/mmax-model-selection/-/blob/gw-distributions/mmms/engine.py#L54) **Update link once merge request has been accepted**.
 
-  * bad population weights? did I have the normalizations messed up?
-    - note that there is excellent agreement when we have a fixed population (all uncertainty is from the EoS), but there is bad agreement when we have uncertainty in the population and in the EoS. This suggests the issue is associated with the marginalization over the population.
-    - add links to where this was/is defined, etc.
+This difference can cause discrepencies between the old and the new code when there is a nontrivial marginalization over the population (i.e., when there is more than a single population sample `\Lambda_p`).
+These can be much larger than the statistical uncertainty from the finite number of Monte Carlo samples (see [below](#comparison-to-previous-results)), but have never been large enough to change the scientific conclusions drawn. Furthermore, [ApJ 915 L5 (2021)](https://iopscience.iop.org/article/10.3847/2041-8213/ac082e) only reported results assuming a fixed population (flat in source-frame component masses) and therefore this issue did not affect the published results.
 
 ## technical note