How to make errors propagation using emcee

Question

I am using emcee in order to analyze some data, related to the SN 1a. Before tackling the real data, I started with simulated data. The results I obtained appear to be good:

Now, I want to calculate the posteriors of some quantities, such as

f(q,j) = j - q^2

using the posteriors I have obtained for q and j. Can anyone help me?

Answer 1

If you have raw posterior samples of j and q, e.g., two numpy arrays j and q, then computing elementwise on those arrays will yield the corresponding samples for the desired variable. In your example, f = j - q**2.

I know options to capture samples of transformed variables of interest directly in the sampler exist in other samplers (Stan, PyMC3). Maybe someone who knows emcee (not me) knows how to do this.

Answer 2

You can use the final samples of the parameters, then get the quantiles of your transformation, for example:

ndim=3 #number of parameters
chains= sampler.chain
samples = chains.reshape((-1, ndim))
q=samples[:,0]
j=samples[:,1]
H=samples[:,2]
f = j - q**2

f now can be represented as a distribution:

import matplotlib.pyplot as plt
import numpy as np
plt.hist(f)
Q=np.quantile(f,q=[0.16,0.5,0.84])

If f gets to a gaussian distribution, then Q are the median Q[1] with the 1 sigma values `Q[0],Q[1]'

You can also use the standard error propagation equation (https://en.wikipedia.org/wiki/Propagation_of_uncertainty):

f=Ej-Eq**2
df = np.sqrt(dj**2 - (2*Eq*dq)**2)

which gives f = 0.9 +/- 0.25 for dq = 0.03, dj = 0.25, Eq=-0.52, Ej=1.18 (from your plot). But due to the high correlations of the parameters this is not the safest approach.