Technical-note-on-implementation.md

-# Introduction
+# LIV parametrization
 
 The bilby implementation follows the implementation of the dispersion implemented in LAL simulation. The parametrization of LIV has the form:
 ```math
@@ -39,7 +39,7 @@ The difference stems from the 2 points:
 1. The correction is parametrized in terms of $` A_{\alpha,eff}`$ instead of $` \lambda_{A,eff}`$
 2. The correction used in lal inference was derived using phase velocity, while this one uses group velocity, following `https://dcc.ligo.org/DocDB/0182/P2200154/001/GW_phase_degeneracies.pdf`.
 
-# Reparametrizing phase velocity correction in terms of A
+## Reparametrizing phase velocity correction in terms of A
 
 1. Route one - plug in the definition of $`\lambda_{\alpha, eff}`$.
 ```math
@@ -62,7 +62,7 @@ Combining equations 28-32 (remembering that they are in units $`G=c=1`$) togethe
 ```
 and then follow the steps of route 1.
 
-# Parametrization in terms of group velocity.
+## Parametrization in terms of group velocity.
 
 Reference `https://dcc.ligo.org/DocDB/0182/P2200154/001/GW_phase_degeneracies.pdf` derives the LIV correction if one uses group velocity instead of phase velocity. The result is eq. 2.9:
 
@@ -74,3 +74,62 @@ The paper uses $`\hbar=1`$ and does not include $'c'$ in the definition of $`D_\
 ```math
 \delta \Phi_{\alpha}(f)= -\frac{\pi D_\alpha [f(1+z)]^{\alpha-1}h^{\alpha-2}}{c}A_{\alpha}= -\frac{\pi D_L h^{\alpha-2}}{c}A_{\alpha,eff}(f )^{\alpha-1}
 ```
+
+# Re-weighting posteriors
+
+Sampling directly on $`A_\alpha`$ would be computationally expensive. The LIV phase term depends on $`D_\alpha`$, which is an integral over redshift. Therefore, if we were to sample on $`A_\alpha`$, we would 1st have to find redshift $`z`$ by numerically inverting $`D_L(z)`$, and then do an integral to find $`D_\alpha(z)`$. Instead, we sample on $`A_{\alpha, eff} =(D_{\alpha}/D_L)(1+z)^{\alpha-1}A_{\alpha}`$, bypassing the problem.
+
+(As a side note we could also try using a lookup table in the likelihood function, with which we could quickly convert between $`D_L`$ and $`D_\alpha`$, but that would require implementing custom bilby likelihood function)
+
+We are interested in the end in how would the posterior look under the prior flat in $`A_\alpha`$ (flat in graviton mass), therefore we need to correct the posterior after the PE run finishes.
+
+## Transforming posteriors
+
+Consider posterior $`p(\theta | d)`$ obtained from sampling in parameters $`\theta`$ and the posterior we would have gotten if we sampled in parameters $`\theta'`$, $`p(\theta' | d)d`$. They are related by the Jacobian:
+
+
+
+```math
+\frac{p(\theta', d)}{p(\theta', d)} = \left| \frac{\partial \theta'_i}{\partial \theta_j} \right| \, . 
+```
+As we transform only one parameter between the 2, the one that parametrizes LIV, the Jacobian matrix reduces to just one component:
+```math
+\frac{p(\theta', d)}{p(\theta', d)} = \left| \frac{\partial \theta'}{\partial \theta} \right| \, , 
+```
+with $`\theta`$ now referring just to the LIV parameter.
+
+For the purpose of reweighting, we don't have to worry about the overall proportionality factor. In the code in this implementation, I account for the following transformations:
+
+1. $`A_{\alpha, eff} \rightarrow A_\alpha`$
+
+```math
+A_\alpha \propto A_{\alpha, eff} \rightarrow \left| \frac{\partial A_\alpha}{\partial A_{ \alpha, eff}} \right| \propto \frac{A_\alpha}{A_{\alpha, eff}}\, , 
+```
+
+2. $`\log(A_{\alpha, eff}) \rightarrow A_\alpha`$
+```math
+\left| \frac{\partial A_\alpha}{\partial \log(A_{\alpha, eff})} \right| \propto \left| \frac{\partial A_\alpha}{\partial \log(A_{\alpha})} \right|  \propto \left| \frac{d A_\alpha}{\frac{1}{A_\alpha} dA_{\alpha}} \right| = |A_\alpha|
+```
+3. $`A_{0, eff} \rightarrow m_g`$
+```math
+A_0 \propto m_g^2 \rightarrow \left| \frac{\partial m_g}{\partial A_{0, eff}} \right| = \left| \frac{\partial m_g}{\partial A_0}\frac{\partial A_0}{\partial A_{0, eff}} \right| \propto \left| \frac{\partial A_0^{\frac{1}{2}}}{\partial A_0}\frac{A_0}{A_{0, eff}} \right| \propto \left| \frac{A_0^{\frac{1}{2}}}{A_{0, eff}} \right|
+```
+This blows up near $`A_0=0`$, so this resampling should better be avoided.
+
+4. $`\log(A_{0, eff}) \rightarrow m_g`$
+```math
+\left| \frac{\partial m_g}{\partial \log(A_{0, eff})} \right| = \left| \frac{\partial m_g}{\partial \log(A_{0})} \right| = \left| \frac{\partial m_g}{\partial A_0}\frac{\partial A_0}{\partial \log(A_{0})} \right| \propto \left| \frac{\partial A_0^{\frac{1}{2}}}{\partial A_0}A_0 \right| \propto \left| A_0^{\frac{1}{2}} \right|
+```
+## Resampling
+
+We can now assign weights $`w`$ to each sample based on the Jacobians computed above, and reweight them using rejection sampling:
+
+1. $`w = w/ w_{max}`$
+2. For each sample, sample point $`x`$ from uniform distribution $`x \in \mathcal{U}_{[0,1]}`$
+3. Keep the sample if $`x<=w`$, refect it otherwise.
+
+## Alternatively, we can use the weights to compute the percentile function for the marginalized distribution,
+```math
+Q=In \left(q, \frac{\sum_{i<j}w}{\sum w}, d_j \right)
+```
+where $`In`$ is interpolation function, $`q`$ is data point to be interpolated and $`d_j`$ are the sorted parameter samples.
\ No newline at end of file