XID+ Example Output Analysis

(This is based on a Jupyter notebook, available in the XID+ package and can be interactively run and edited)

This notebook provides some example code for basic analysis of the XID+ outputs, including:

  1. Loading up output
  2. Creating Posterior replicated maps and animations
  3. Creating marginalised posterior plots
  4. Creating Bayesian p-value maps

Import required modules


In [1]:
import pylab as plt
%matplotlib inline


import numpy as np
import xidplus
from xidplus import moc_routines
output_folder='./'

Load up posterior output from XID+


In [2]:
priors,posterior=xidplus.load('test.pkl')

In order to compare how good our fit is, its often useful to look at original map. There is a routine within XID+ that makes the original fits map from the data stored within the prior class. Lets use that to make the SPIRE maps for the region we have fit.

Now lets use the Seaborn plotting package and APLpy package to view those maps, plotting the sources we have fit on top of those maps.


In [3]:
figs,fig=xidplus.plot_map(priors)


Posterior replicated data

We can use each sample we have from the posterior, and use it to make a replicated map, including simulating the instrumental noise, and the estimated confusion noise. You can think of these maps as all the possible maps that are allowed by the data.

NOTE: You will require the FFmpeg library installed to run the movie


In [4]:
xidplus.replicated_map_movie(priors,posterior,50)


Out[4]:

Joint Posterior Analysis of sources


In [5]:
#Select source you want to plot joint distribution
s1=2
s2=18

Sources 2 and 18 are close together, lets look at their joint posterior probabiity distribution. The code below is an example of how you can plot the joint posterior probability density function for the $250 \mathrm{\mu m}$ flux, and an inset of the real map.


In [6]:
import aplpy
import seaborn as sns
sns.set(color_codes=True)
import pandas as pd
sns.set_style("white")
import xidplus.posterior_maps as postmaps



labels=[r'Source 1 $250\mathrm{\mu m}$ flux  (mJy)',r'Source 2 $250\mathrm{\mu m}$ flux (mJy)']
df = pd.DataFrame(posterior.samples['src_f'][:,0,[s1,s2]],columns=labels)
g = sns.PairGrid(df,size=5)
g.map_diag(sns.kdeplot,c='Red')
g.map_lower(sns.kdeplot, cmap="Reds",alpha=0.8,n_levels=10,normed=True, shade=True,shade_lowest=False)
g.set(ylim=(0,40))
g.set(xlim=(0,40))

g.axes[0,1].spines['bottom'].set_color('white')
g.axes[0,1].spines['left'].set_color('white')
cmap=sns.cubehelix_palette(8, start=.5, rot=-.75,as_cmap=True)

real_250 = aplpy.FITSFigure(postmaps.make_fits_image(priors[0],priors[0].sim)[1],figure=g.fig,subplot=(2,2,2))
real_250.show_colorscale(cmap=cmap)

real_250.show_markers(priors[0].sra, priors[0].sdec, edgecolor='black', facecolor='black',
                marker='o', s=40, alpha=0.5)
real_250.recenter(priors[0].sra[s1], priors[0].sdec[s1], radius=0.01)

real_250.add_label(priors[0].sra[s1], priors[0].sdec[s1]+0.0005, 1, relative=False,size=20,color='white')
real_250.add_label(priors[0].sra[s2], priors[0].sdec[s2]-0.0010, 2, relative=False,size=20,color='white')
real_250.tick_labels.set_xformat('dd.dd')
real_250.tick_labels.set_yformat('dd.dd')


real_250.add_colorbar(axis_label_text=r'$250\mathrm{\mu m}$ flux  (mJy)')


INFO: Auto-setting vmin to -1.522e+01 [aplpy.core]
INFO: Auto-setting vmax to  3.332e+01 [aplpy.core]

Posterior Predictive checking and Bayesian P-value maps

When examining goodness of fits, the typical method is to look at the residuals. i.e. $\frac{data - model}{\sigma}$. Because we have distribution of $y^{rep}$, we can do this in a more probabilisitic way using posterior predictive checks. For more information on posterior predictive checks, Gelman et al. 1996 is a good starting point.

For our case, the best way to carry out posterior predictive checks is to think about one pixel. We can look at where the real flux value for our pixel is in relation to the distribution from $y^{rep}$.


In [12]:
from xidplus import posterior_maps as postmaps

rep_maps=postmaps.replicated_maps(priors,posterior)

import matplotlib as mpl
sns.set_style("white")

fig=plt.figure(figsize=(5,5))


# This is  the colormap I'd like to use.
cm = sns.diverging_palette(220, 20, as_cmap=True)

# Get the histogramp
Y,X = np.histogram(rep_maps[0][20,:], 25, normed=1)
#C = [cm(((x-X.min())/x_span)) for x in X]
C = [cm(((((x-np.mean(rep_maps[0][20,:]))/np.std(rep_maps[0][20,:]))+6)/12.0)) for x in X]


plt.bar(X[:-1],Y,color=C,width=X[1]-X[0])
plt.xlabel('Pixel Flux mJy')
plt.ylabel('p(pixel flux)')
plt.axvline(5, linestyle='--')
plt.axvline(-10,linestyle=':')
plt.annotate('flux in map that \n model cannot explain',xy=(5, 0.01),  xycoords='data',
            xytext=(5, 0.1), textcoords='data',rotation='vertical',size='small')
plt.annotate('too much flux in model \n compared to map',xy=(-10, 0.01),  xycoords='data',
            xytext=(-10, 0.1), textcoords='data',rotation='vertical',size='small')

#ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
ax1 = fig.add_axes([0.95, 0.1, 0.05, 0.8])

norm = mpl.colors.Normalize(vmin=-6, vmax=6)
cb1 = mpl.colorbar.ColorbarBase(ax1, cmap=cm,
                                norm=norm,
                                orientation='vertical')
cb1.set_label('$\sigma$')


We can calculate fraction of $y^{rep}$ samples above and below real map value. This is often referred to as the Bayesian p-value and is telling us the probability of drawing the real pixel value, from our model which has been inferred on the data. This is tells us if the model is inconsistent with the data, given the uncertianties in parameters and data.

  • $\sim 0.5$ means our model is consistent with the data
  • 0.99 or 0.01 means model is missing something.

We can convert this to a typical '$\sigma$' level and create map versions of these Bayesian p-values:


In [13]:
figs, fig=xidplus.plot_Bayes_pval_map(priors, posterior)


Red indicates the flux value in the real map is higher than our model thinks is possible. This could be indicating there is a source there that is not in our model. Blue indicates the flux in the real map is lower than in our model. This is either indicating a very low density region or that too much flux has been assigned to one of the sources.

Creating Catalogues

We can also create catalogues from the posterior probability density function


In [14]:
import xidplus.catalogue as cat

In [16]:
SPIRE_cat=cat.create_SPIRE_cat(posterior,priors[0],priors[1],priors[2])
SPIRE_cat.writeto('test.fits',overwrite=True)

Check table

Lets read in table with Astropy table


In [20]:
from astropy.table import Table

In [21]:
catalogue=Table.read('test.fits')

In [22]:
catalogue


Out[22]:
<Table length=51>
HELP_IDRADecF_SPIRE_250FErr_SPIRE_250_uFErr_SPIRE_250_lF_SPIRE_350FErr_SPIRE_350_uFErr_SPIRE_350_lF_SPIRE_500FErr_SPIRE_500_uFErr_SPIRE_500_lBkg_SPIRE_250Bkg_SPIRE_350Bkg_SPIRE_500Sig_conf_SPIRE_250Sig_conf_SPIRE_350Sig_conf_SPIRE_500Rhat_SPIRE_250Rhat_SPIRE_350Rhat_SPIRE_500n_eff_SPIRE_250n_eff_SPIRE_500n_eff_SPIRE_350Pval_res_250Pval_res_350Pval_res_500
degreesdegreesmJymJymJymJymJymJymJymJymJymJy/BeammJy/BeammJy/BeammJy/BeammJy/BeammJy/Beam
str27float64float64float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32
1891150.747114622.019372291493.74847.663191.181642.143744.784480.6541443.956588.780071.13176-2.55132-4.50904-6.236811.575671.541542.358091.001250.9996510.9996942000.02000.02000.00.00.0010.0
1896150.746556442.013817344354.600056.36972.842673.658665.486891.850561.928264.578170.519253-2.55132-4.50904-6.236811.575671.541542.358091.001530.9992811.000562000.02000.02000.00.00.0030.0
2640150.7426721672.0293733060413.998418.88439.1746312.898115.48719.418061.869514.19930.517889-2.55132-4.50904-6.236811.575671.541542.358091.001921.001161.000282000.02000.02000.00.00.0080.0
5128150.7526802822.036591142423.179355.001121.515381.175022.57950.3161554.843748.411981.72842-2.55132-4.50904-6.236811.575671.541542.358090.9992810.9997250.9990162000.02000.02000.00.00.0170.001
5129150.7410062882.03326257014.140865.904872.338281.572663.093270.5151442.126924.432630.572698-2.55132-4.50904-6.236811.575671.541542.358090.9996270.9993960.9994932000.02000.02000.00.00.0180.0
5130150.7471233192.0399264171415.16316.773413.559918.350819.890316.837814.38617.360211.346-2.55132-4.50904-6.236811.575671.541542.358090.9989790.9985970.9990772000.02000.02000.00.00.0030.001
5159150.7543466032.0338129052415.013916.807713.05777.175889.162444.820384.098727.768031.15535-2.55132-4.50904-6.236811.575671.541542.358090.999181.001120.9993272000.02000.02000.00.0020.0040.001
6067150.7554566052.029923855242.47554.139461.05965.077327.130543.021582.4075.264340.697328-2.55132-4.50904-6.236811.575671.541542.358091.001530.9998710.9992672000.02000.02000.00.0020.0230.002
6728150.7315579422.035488248892.275464.036990.8145941.282562.656320.3771021.257022.829280.311752-2.55132-4.50904-6.236811.575671.541542.358090.9991940.9990040.9992942000.02000.02000.00.0010.00.002
.................................................................................
54557150.7471160312.022705389674.166086.520532.010111.809363.885620.4812753.853728.231031.19911-2.55132-4.50904-6.236811.575671.541542.358090.9991270.9989050.9995052000.02000.02000.00.00.0030.0
56696150.7215561732.043824810930.884122.16440.2407980.9493732.360940.2554321.200043.066910.313454-2.55132-4.50904-6.236811.575671.541542.358090.9993510.9989611.001712000.02000.02000.00.0140.0020.008
56700150.7293370912.041599802171.482472.863670.4490761.045822.223370.3227391.391983.187430.412141-2.55132-4.50904-6.236811.575671.541542.358090.9999420.9998220.998372000.02000.02000.00.0010.00.004
57200150.7404581262.0515948904110.47412.7447.88944.536386.913962.202412.684145.658350.824703-2.55132-4.50904-6.236811.575671.541542.358090.9986460.998650.9990782000.02000.02000.00.00.0010.005
58792150.7393341292.022152857781.535842.89060.5367871.040912.183790.3118711.339512.953360.337598-2.55132-4.50904-6.236811.575671.541542.358090.9994911.000850.9995122000.02000.02000.00.00.0020.0
62679150.7510158472.043813518312.448074.682180.7555752.667075.367110.9643113.189247.197390.878154-2.55132-4.50904-6.236811.575671.541542.358090.9998910.9996941.000342000.02000.02000.00.00.00.004
62696150.7165472932.028272101750.9001131.877450.2733191.051672.089770.3232882.65564.829450.958473-2.55132-4.50904-6.236811.575671.541542.358090.9986411.000140.9988872000.02000.02000.00.0010.0680.226
64788150.7571178952.015479799417.40538.999435.699526.234988.095414.302851.509783.650690.413473-2.55132-4.50904-6.236811.575671.541542.358090.9991240.998980.9992342000.02000.02000.00.00.0030.004
64789150.7226685432.045490955680.6934141.731350.1719591.009552.435270.2540881.504513.613530.414402-2.55132-4.50904-6.236811.575671.541542.358090.9995811.000641.000242000.02000.02000.00.3190.0010.002
64790150.7460060082.026594435443.666455.708631.700929.7518511.8367.742272.950045.881840.874799-2.55132-4.50904-6.236811.575671.541542.358091.000891.000120.9996322000.02000.02000.00.00.0040.0