XD is a code for estimating an underlying probability density of some physical features in the presence of (understood!) experimental covariances. In short, given data with known noise properties, XD "deconvolves" the known covariances from the noisy distribution, producing an estimate of data that would have been obtained with an "ideal" apparatus.
In the lQSO classification problem, we want to estimate the following using Bayes' rule:
\begin{align} P(lQSO|features) = \frac{P(features|lQSO)P(lQSO)}{P(features|lQSO)P(lQSO) + P(features| \neg lQSO)P( \neg lQSO)} \end{align}$P(lQSO)$ is a prior probability of an object in our catalog being a lensed quasar.
The likelihood function, $P(features|lQSO)$, depends on both underlying physics and the noise properties of our experiment. Beacuse many of the lQSOs we expect from PS1 will have noisy measurements, and because this noise may lead to fake detections, we need to use XD to evaluate this likelihood function in terms of underlying physics only.
Unfortunately, XD is fundamentally a gaussian mixture model (GMM), and is therefore most efficient at modeling distributions which are close to multivariate gaussians. This means that number counts of astronomical objects (which obey a power law scaling), though physically important, will require the XD model to have many (>20) gaussian components, which is highly inefficient. We have two options for dealing with non-gaussian features:
The second is the choice made by Bovy et al in the SDSS quasar targeting paper. In each bin, the likelihood function is now represented as:
\begin{align} P(features|lQSO) = P(\text{gaussian features }|\text{ i-band mag, lQSO}) P(\text{i-band mag }|lQSO) \end{align}where the first term can now be estimated using XD, and the second term is another "prior", estimated from:
Similar reasoning applies to the second denominator term, where $P( \neg lQSO)$ is just $1-P(lQSO)$ and $P(features| \neg lQSO)$ is estimated:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
dataLoc = '../../'
labels = ['ID','mag_model_i','g-r', 'r-i', 'i-z', 'WISE1', 'WISE2' ]
pQSO = np.loadtxt(dataLoc+'pQSO/pSDSScolmag.txt')
lQSO = np.loadtxt(dataLoc+'lQSO/SDSScolmag.txt')
sinQSO = np.loadtxt(dataLoc+'sinQSO/sSDSScolmag.txt')
unlQSO = np.loadtxt(dataLoc+'unlQSO/nlSDSScolmag.txt')
unlQSO[:,5:7] = -unlQSO[:,5:7] #bug in WISE magnitudes for this file
duds = np.concatenate((pQSO,unlQSO,sinQSO),axis=0)
data = np.concatenate((lQSO,duds),axis=0)
print lQSO.shape, duds.shape, data.shape
numPts = data.shape[0]
data = data[:,2:] #don't use IDs or i-band magnitudes
In [6]:
from astroML.density_estimation import XDGMM
from astroML.plotting.tools import draw_ellipse
import triangle
nFeatures = data.shape[1]
#generate random covariance matrix; perturb data
perturbingSigma = .2
z = zip(np.random.normal(scale=perturbingSigma,size=nFeatures),np.random.normal(scale=perturbingSigma,size=nFeatures))
covmat = np.cov(z)
perturbations = np.random.multivariate_normal(np.zeros(nFeatures),covmat,size=data.shape[0])
noisyData = data + perturbations
#XD function needs to be delivered a list of covariance matrices
covmatlist = np.dstack((covmat,covmat))
while covmatlist.shape[2] < data.shape[0]:
covmatlist = np.dstack((covmatlist,covmat))
# fit the model
xdgmm = XDGMM(n_components=7)
xdgmm.fit(noisyData, covmatlist.transpose())
logp=xdgmm.logprob_a(noisyData, covmatlist.transpose()) # evaluate probability
resampledXD = xdgmm.sample(data.shape[0])
In [25]:
_ =plt.hist(np.sum(np.exp(logp),axis=1),bins=20)
print data.shape
print np.max(np.exp(logp))
In [10]:
#plot the resampled data drawn from the XD density model
text = ['g-r', 'r-i', 'i-z', 'WISE1', 'WISE2' ]
fig1 = triangle.corner(data,labels=text)
fig1.suptitle('(Noisy) Data plotted over \'Truth\'')
_ = triangle.corner(noisyData, fig=fig1,color='r')
fig2 = triangle.corner(noisyData,color='r',labels=text)
fig2.suptitle('Corrected Data plotted over Noisy data')
_ = triangle.corner(resampledXD,fig=fig2,color='b')
fig3 = triangle.corner(data,labels=text)
fig3.suptitle('Corrected Data plotted over \'Truth\'')
_ = triangle.corner(resampledXD,fig=fig3,color='b')
fig4 = triangle.corner(data,labels=text)
fig4.suptitle('Everything')
_ = triangle.corner(noisyData, fig=fig4,color='r')
_ = triangle.corner(resampledXD,fig=fig4,color='b')
In [26]:
"""
class Dataset(object):
__init__(self)
#fits header properties
plot(self):
#make triangle plot; overlay on current triangle plot, if any.
return
"""
# import PS1QLS
from astroML.density_estimation import XDGMM
from astroML.plotting.tools import draw_ellipse
import triangle
class Classifier(object):
def __init__(self,algorithm='XD'):
if algorithm == 'XD':
self.xdmodel = XDGMM(n_components=7)
else:
print 'Not implemented yet!'
return
def train(self,training,covmatlist):
self.xdmodel.fit(training, covmatlist.transpose())
return
def test(self,test):
self.xdprobs=np.sum(np.exp(self.xdmodel.logprob_a(test, covmatlist.transpose())),axis=1)
return
def report(self,truth):
return
RF = Classifier(algorithm='XD')
RF.train(noisyData,covmatlist) #
RF.test(noisyData) # calculates classification probabilities
plt.hist(RF.xdprobs,bins=20)
#RF.report() # make plots, report statistics etc.
#RF.save() # pickles itself
#RF.score(testtruth) # plots an ROC curve
Out[26]:
In [28]:
print RF.xdprobs[100:200]
In [2]:
from astroML.density_estimation import XDGMM
X = np.random.normal(size=(1000, 1)) # 1000 pts in
ax1 = plt.subplot(211)
plt.title('XD Before and After')
_ = plt.hist(X,bins=20)
Xerr = np.random.random((1000, 1, 1)) # 1000 1x1 covariance matrices
xdgmm = XDGMM(n_components=2)
xdgmm.fit(X, Xerr) # fit the model
logp=xdgmm.logprob_a(X, Xerr) # evaluate probability
X_new = xdgmm.sample(1000)
plt.subplot(212,sharex=ax1)
_ = plt.hist(X_new,bins=20)
print 'Sigma before = ', np.std(X), 'Sigma after = ', np.std(X_new)
In [4]:
print np.exp(xdgmm.logprob_a(np.zeros((1,1)),.1*np.ones((1,1,1))))
In [3]:
import brewer2mpl
from IPython.display import display
bmap = brewer2mpl.get_map('Set1', 'qualitative', 4)
col = bmap.mpl_colors
def compare(x,y):
plt.figure()
plt.scatter(pQSO[:,x],pQSO[:,y],c=col[1],alpha=.75,label='QSO+QSO Alignment')
plt.scatter(lQSO[:,x],lQSO[:,y],c=col[0],alpha=.75,label='Lensed QSO')
plt.scatter(sinQSO[:,x],sinQSO[:,y],c=col[2],alpha=.75,label='Single QSO')
plt.scatter(unlQSO[:,x],unlQSO[:,y],c=col[3],alpha=.75,label='QSO+LRG Alignment')
plt.xlabel(labels[x])
plt.ylabel(labels[y])
plt.legend(numpoints=1,bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
In [13]:
interact(compare,x=(0,6),y=(0,6))