Model

We would like to infer the distance to twelve velocity slices of a spectral cube, covering the molecular cloud complex Cepheus. We need two datasets. First, we need the photometry of stars towards Cepheus, which have been used in Green et al. 2015 to derive joint posteriors on distance and reddening for each star. Second, we need the intensity of carbon monoxide emission in each velocity slice of a spectral cube towards this same region; we use the Galaxy-wide CO survey from Dame et al. 2001.

This is what we ultimately wish to infer:

$$p(\alpha \mid m, I) $$

where $$\alpha=(d_1, d_2, ... , d_{12}; c_1,c_2,...,c_{12})$$ and...
$d$=distance to slice (12 d's for 12 velocity slices) $c$=gas to dust coefficient for slice (12 c's for 12 velocity slices) $I$= gas emission intensities that will act as dust templates, and which will be converted to extinction with a gas-to-dust coefficient
$m$= set of stellar photometry for all stars

From Bayes Rule, we know that, up to a normalizing constant:

$$p(\alpha \mid m, I) \propto p(m, I \mid \alpha) \times p(\alpha) $$

with the first term being our likelihood function and the second terms being our prior function; we're assuming a flat prior for $\alpha$.

As explicated in Green et al. 2015, we can show that our likelihood function for each star is simply the line integral over the joint posterior on distance and reddening for that star, taken along some specified reddening profile (cumulative reddening as a function of distance). We introduce a nuisance parameter, $\theta$, the stellar type, and we re-express the likelihood function as an integral over the distance modulus $\mu$ and the stellar type $\theta$:

$$p(m, I \mid \alpha) \propto \int d\mu \; d\theta \; p(m, I, \mu, \theta \mid \alpha)$$

We can factor this into two terms using the chain rule for probability:

$$ \propto \int d\mu \; d\theta \; p(m, I \mid \mu, \theta, \alpha) p(\mu, \theta \mid \alpha)$$

And, if we assume that $\mu$ and $\theta$ are both independent of the reddening profile, defined by the $\alpha$ coefficients, we get

$$\propto \int d\mu \; d\theta \; p(m, I \mid \mu, \theta, E(\mu;\alpha)) \; p(\mu, \theta)$$

The integrand in the last equation is basically (up to a normalizing constant) a prior times a likelihood, so our likelihood for one star is essentially the integral of the posterior density $p(\mu,E, \theta \mid m)$. We marginalize out stellar type, so we're left with an the integral over $\mu$. The posterior density for each star has already been precomputed in Green et al. 2015. Given that the likelihood for all stars is given by the product of the individual likelihoods, our posterior for $\alpha$ is given by:

$$p(\alpha \mid I, m) \propto p(\alpha) \times \prod_j \int p(\mu_j, E_j(\mu_j; \alpha) \mid m_j, I) \; d\mu_j$$

Now we need a model for the reddening profile, which defines the path of our line integral through the stellar posterior arrays. We parameterize the line-of-sight reddening profile towards an individual star, $E_j$, with the following model:

$$E_{j} (\mu, \alpha) = \sum_{k \rightarrow d_k < \mu} (I_{j,k})\times c_k$$

with $I_{j,k}$= intensity of gas emission in the pixel corresponding to the jth star in the kth velocity slice and $c$= gas to dust conversion coefficient for the kth velocity slice. We're summing over the emission in all the velocity slices such that the distance to the velocity slice is less than the distance to our star $\mu$. The form of the reddening profile will be piecewise linear. Our coefficient converts CO intensity in units of Kelvin to reddening in units of magnitudes. While we do let the coefficient float, we need to set the walkers off at a reasonable value from the literature. So, to derive this, we convert from CO intensity to atomic hydrogen column density (Dame et al. 2001) and then from atomic hydrogen column density to reddening (Schlegel, Finkbeiner, and Davis et al. 1998). Finally, we multiply by the width of each spectral channel in $\rm km \; s^-1$ to get units of mag/K:

$$\rm c=(3.6\times 10^{20} \; H \; cm^{-2} \; K^{-1} \; km^{-1} \; s) \times (\frac{1 \; mag \; E(B-V)}{8 \times 10^{21} \; H \; cm^{-2}}) \times (\frac{\rm 1.3 \; km\; s^{-1} }{channel})$$$$ = \frac{\rm 0.06 \; mag}{K}$$

Now that we know the posterior we're sampling from, along with a model for the reddening profile, let's check out how our code works on simulated data...

Tests on Simulated Data

Now we need to do some code testing using simulated data, to ensure that our code is functioning properly. To do this, we need to produce simulated data such that we know what the "true" answer is for the parameter estimation problem we are trying to solve. To do this, we are going to create a reddening "profile" or a plot of the cumulative reddening as a function of distance. In this profile, there are going to be two clouds along the line of sight, corresponding to two jumps in reddening in our reddening profile. We are going to have our reddening profile look like this:

From the reddening profile, we can tell that this is a two cloud model towards a specific line of sight. We see that there are no clouds (and thus no reddening) between a distance modulus of 4 and 7.75. Then, we see that there's a cloud of dust at 7.75 which adds 2 magnitudes of reddening to our profile. Then we see that there's a second dust cloud at a distance modulus of 14 that adds 3 magnitudes of reddening to our profile, bringing our cumulative reddening at distance modulus 14 and beyond up to 5 magnitudes.

So, now we need to produce some stellar posterior data (700x120 arrays of reddening vs. distance) to match our contrived profile. To do this, we are going to randomly sprinkle 50 stars along the line of sight. Since we are working with stellar posterior arrays, we are going to mentally convert everything from actual distance and reddening to distance and reddening bins in the stellar posterior arrays. Axis 0 of the array (the vertical axis) has a reddening range from 0-7 magnitudes, $\Delta_{reddening}$=0.01 magnitudes (total of 700 bins). Axis 0 of the array (the horizontal axis) has a distance modulus range from 4-19, so $\Delta_{distmod}$=0.125 (total of 120 bins). Our clouds in our contrived reddening profile correspond to distance modulus bins of 30 and 80, and reddening bins of 200 and 500. So now we randomly sprinkle 50 stars along the distance axis of the array, and then look up the reddening those stars should have based on our reddening profile. To produce the posterior for a single star we draw 5000 random samples from bivariate normal distribution with a mean centered at the randomly generated distance and the corresponding reddening value given by the reddening profile. We let the standard deviation of the bivariate normal distribution be 0.1 in both distance and reddenings, resulting in "puffy" clouds. We then use the samples to populate a 2D histogram, which is now our 700x120 stellar posterior array for a single star. We repeat this 49 more times for 49 more stars. This is how we code that:

We've plotted below what the posterior would look like if you stacked all the individual stellar posteriors on top of eachother (i.e. you summed along the depth axis). See that it looks similar to our contrived reddening profile, which is what we want! Note that the x and y axes correspond to the bin indices. When we integrate over distance modulus in our likelihood function, we only care about the reddening "ledges". We'll perform this integral from distance modulus=4 to distance modulus=19, over the entire length of the stellar posterior array for each individual star.

Now we are going to make some fake CO (carbon monoxide) data to go with our fake stellar posteriors. We have twelve CO intensity values for each star (for twelve CO velocity slices). The set of intensity values will be the same for each star, since they are all presumed to be along the same line of sight. Our fake reddening profile corresponds to two clouds (two CO emission features) at two different velocity slices. The first cloud has a reddening of 2 magnitudes, producing a cumulative reddening of 2 magnitudes. The second cloud has a reddening of 3 magnitudes, producing a cumulative reddening of 5 magnitudes (3 from itself and 2 from the other cloud). Assuming gas to dust coefficients to one, this means that we need to set one slice to CO intensity ("CO_I") value of 2 and the second slice to a CO intensity value of 3. All the other velocity slice intensities will be set to zero, because we are assuming there are only two clouds in all the slices.

It doesn't really matter which two slices we put our clouds in (it only matters what distances and coefficients our parameter estimation problem gives us for the two slices). Thus, let's put the first cloud (CO_I=2) in velocity slice 7 and our second cloud (CO_I=3) in velocity slice 11. So for a single star, our CO intensities for velocity slices 1-12 look like this: [0,0,0,0,0,0,2,0,0,0,0,3,0]. Now we have a stellar posterior array and a CO intensity array for every star. We write all this information out into a "simulated_data.h5" file:

Based on our stellar posterior and CO data, we know what the values of d7, d11,c7, and c11 SHOULD be. D7 needs to be before D11, at a distance modulus of 7.75 and D11 needs to be at a distance modulus of 14. Both c7 and c11 should have values near 1, because that's what we assumed when created our reddening profile and assigning CO intensity values. Having those velocity slices correspond to those distances and coefficients will produce the highest value of the line integral along our reddening profile (see initial reddening profile figure), which is essentially our likelihood function. The distances and coefficients to the other ten slices cannot be constrained because their intensities are set to zero and thus contribute nothing to the reddening profile.

Because our distributions are likely to be multi-modal, we are going to use Parallel-Tempering MCMC. To ensure that we aren't getting the "right" answer by chance, we are going to randomly set our walkers off within a reasonable distance range surrounding each cloud. For the d7 parameter, we are going to set our walker off somewhere between a distance modulus of 8 and 11, for our d11 parameter we are going to set our walker off somewhere between a distance modulus of 14 and 17. This way, we are sure not to initialize our walkers at an area of negative infinity probability, while still making it somewhat difficult for them the find their way back. Our final answers don't have to be exact (because we'd need a lot more simulated stars along the line of sight to make it perfect) but it does have to be close. Moreover, our distance estimates should get better if we increase the number of stars and should get worse if we increase the standard deviation of the bivariate Gaussian's we used to create the probability surfaces in our stellar posterior array. Let's see what happens...

Results

Success! We see that we returned a value of for d7, d11, c7, and c11 are close to the true values given by 7.75, 14.00, 1.0, and 1.0. Thus, our PTSampler produced accurate results, well within the estimated uncertainties.

Please note: if you're rerunning this notebook, these values will vary from trial to trial, and will depend on how many stars you simulate, along with the random distance generated along the sight for each star; every time you run the entire notebook and recreate the simulated_data.h5 file, you will get a slightly different reddening profile with stars at different distances along the $\mu$ axis of the posterior array.

Finally, let's check out our reddening profile, using the 50th quantile samples returned by the MCMC chains. To see how well we did, let's overplot our reddening profile on top of our stacked stellar posterior surface.

Looks good! Let's move on to real data in tutorial.ipynb