Selection of Burst-like Transients and Stochastic Variables Using Multi-Band Image Differencing in the Pan-STARRS1 Medium-Deep Survey

1.2 MB·English

by  S. Kumar

Draft version January 8, 2015
PreprinttypesetusingLATEXstyleemulateapjv.08/13/06
SELECTION OF BURST-LIKE TRANSIENTS AND STOCHASTIC VARIABLES USING MULTI-BAND IMAGE
DIFFERENCING IN THE PAN-STARRS1 MEDIUM-DEEP SURVEY
S. Kumar1, S. Gezari1, S. Heinis1, R. Chornock2, E. Berger2, A. Rest3, M.E. Huber4, R.J. Foley6,10, G.
Narayan4, G.H. Marion4, D. Scolnic7, A. Soderberg2, A. Lawrence5, C.W. Stubbs2, R.P. Kirshner2, A.G. Riess7,
S.J. Smartt8, K. Smith8, W.M Wood-Vasey9, W. S. Burgett4, K. C. Chambers4, H. Flewelling4, N. Kaiser4, N.
Metcalfe3, P. A. Price4, J. L. Tonry4, & R. J. Wainscoat4
Draft version January 8, 2015
5
1 ABSTRACT
0
We present a novel method for the light-curve characterization of Pan-STARRS1 Medium Deep
2
Survey(PS1MDS)extragalacticsourcesintostochasticvariables(SV)andburst-like(BL)transients,

n usingmulti-bandimage-differencingtime-seriesdata. Weselectdetectionsindifferenceimagesassoci-
a atedwithgalaxyhostsusingastar/galaxycatalogextractedfromthedeepPS1MDSstackedimages,
J
and adopt a maximum a posteriori formulation to model their difference-flux time-series in four Pan-

6 STARRS1 photometric bands g ,r ,i , and z . We use three deterministic light-curve models
P1 P1 P1 P1

tofitburst-liketransients;aGaussian,aGammadistribution,andananalyticsupernova(SN)model,
]
M and one stochastic light curve model, the Ornstein-Uhlenbeck process, in order to fit variability that
is characteristic of active galactic nuclei (AGN). We assess the quality of fit of the models band-wise
I source-wise, using their estimated leave-out-one cross-validation likelihoods and corrected Akaike in-
.
h formation criteria. We then apply a K-means clustering algorithm on these statistics, to determine
p thesourceclassificationineachband. Thefinalsourceclassificationisderivedasacombinationofthe
- individual filter classifications, resulting in two measures of classification quality, from the averages
o
across the photometric filters of 1) the classifications determined from the closest K-means cluster
r
t centers, and 2)the square distancesfrom theclustering centers inthe K-means clusteringspaces. For
s
a verification set of active galactic nuclei and supernovae, we show that SV and BL occupy distinct
a
[ regions in the plane constituted by these measures. We use our clustering method to characterize
4361 extragalactic image difference detected sources in the first 2.5 years of the PS1 MDS, into 1529
1
BL, and 2262 SV, with a purity of 95.00% for AGN, and 90.97% for SN based on our verification
v
sets. Wecombineourlight-curveclassificationswiththeirnuclearoroff-nuclearhostgalaxyoffsets,to
4
definearobustphotometricsampleof1233activegalacticnucleiand812supernovae. Withthesetwo
1
samples, we characterize their variability and host galaxy properties, and identify simple photometric
3
priors that would enable their real-time identification in future wide-field synoptic surveys.
1
0 Subject headings: methods: statistical, techniques: photometric, (galaxies:) quasars: general
.
1
0
1. INTRODUCTION well as in real-time data. Recent increases in computa-
5
tional resource availability and efficiency have enabled
1 With the advent of the LSST era, human-intervened
near-complete automated transient discovery in large
: classification (with the exception of citizen science) will
v surveys (Bloom et al. 2012). Machine-learning methods
become untenable, requiring automated source identifi-
i areslowlyreplacinghumanjudgementfortransientclas-
X cation in large volumes of data in archival catalogs, as
sification in real-time, as well as in large survey catalogs
ar 1Department of Astronomy, University of Maryland, Stadium (Richards et al. 2011; Sanders et al. 2014; Pichara &
Drive,CollegePark,MD21224,USA Protopapas 2013). The knowledge of prior event types
2Harvard-Smithsonian Center for Astrophysics, 60 Garden makes it possible to look for specific events in the data
Street,Cambridge,Massachusetts02138,USA
with a high degree of completeness and efficiency us-
3Space Telescope Science Institute, 3700 San Martin Drive,
ing time-variability (Bailer-Jones 2012; Butler & Bloom
Baltimore,Maryland21218,USA
4InstituteforAstronomy,UniversityofHawaii,2680Woodlawn 2011; Schmidt et al. 2010; Choi et al. 2013), color based
Drive,Honolulu,Hawaii96822,USA selection (Richards et al. 2009), multi-wavelength cata-
5Institute for Astronomy, University of Edinburgh Scottish
log associations (Mendez et al. 2013), and host-galaxy
Universities Physics Alliance, Royal Observatory, Blackford Hill,
EdinburghEH93HJ,UK properties (Foley & Mandel 2013). Also, generalized
6Astronomy Department, University of Illinois at Urbana- automated machine classification algorithms based on
Champaign,1002WestGreenStreet,Urbana,IL61801,USA random-forest methods (Lo et al. 2014), support vec-
7Department of Physics and Astronomy, Johns Hopkins Uni-
tor machines and naive Bayes estimates (Mahabal et al.
versity, 3400 North Charles Street, Baltimore, Maryland 21218,
USA 2008), and sparse matrix methods (Wozniak et al. 2013)
8Astrophysics Research Centre, School of Mathematics and that use a number of photometric and non-photometric
Physics,Queen’sUniversityBelfast,BelfastBT71NN,UK features have been demonstrated to achieve classifica-
9Pittsburgh Particle Physics, Astrophysics, and Cosmology
tions with very high purity.
Center, Department of Physics and Astronomy, University of
Pittsburgh, 3941 O’Hara Street, Pittsburgh, Pennsylvania 15260, As the number of detected transients grows very large
USA in wide-field time domain surveys, complete spectro-
10Department of Physics,University of Illinois Urbana- scopic follow up becomes impossible due to limited re-
Champaign,1110W.GreenStreet,Urbana,IL61801USA
2
sources and faint magnitude limits. Classification meth-
ods using time-series data alone are favorable, and have TABLE 1
Pan-STARRS1 Medium Deep
been applied in the past to a broad range of sources;
Survey Field Centers
Choi et al. (2013) discuss the identification of AGN
via damped-random walk parameterization of image-
··· J2000 J2000
differencing light curves, Andrae et al. (2013) on the Field RA Declination
applicability of single and multiple Ornstein-Uhlenbeck ··· Degrees Degrees
(OU)processestoAGNlightcurves,andButler&Bloom
MD01 02h23m30s −04deg15(cid:48)
(2011) on the separation of AGN from variable stars in MD02 03h32m24s −27deg48(cid:48)
photometric surveys through damped-random walk pa- MD03 08h42m22s 44deg19(cid:48)
rameterization. For supernovae (SNe), Kessler et al. MD04 10h00m00s 02deg12(cid:48)
(2010) discusses various photometric template fitting MD05 10h47m40s 58deg04(cid:48)
methods that enable their identification with particular MD06 12h20m30s 47deg07(cid:48)
SN classes. MD07 14h14m48s 53deg04(cid:48)
MD08 16h11m08s 54deg57(cid:48)
Bailer-Jones (2012) gives a powerful method to pick
MD09 22h16m45s 00deg16(cid:48)
suitablemodelsfordeterministicorstochasticlightcurves MD10 23h29m14s 00deg25(cid:48)
using leave-one-out cross-validation. This method can
be used to determine source type, based on the likeli- properties of extragalactic variables and transients, and
hoods of the various models for a given lightcurve. So combine our light-curve classifications with host galaxy
far, the applicability of selecting sources based on time- offsets in order to define a robust photometric sample
series modeling has been limited to single-band detec- of SNe and AGN; and in §6 we describe the source and
tions (Choi et al. 2013), or have typically used magni- host galaxy properties of our variability/offset selected
tude time-series data (Butler & Bloom 2011), which are SN and AGN samples, and report on priors that can aid
undefined for negative difference-fluxes. Also, computa- in their real-time classification in future surveys.
tional limitations typically have lead to the use of only
singlemodelsaspredictorsforclass,oronlyusingsimple 2. THEPAN-STARRS1SURVEYANDTRANSIENTALERT
DATABASE
statistical criteria for model assessment, rendering clas-
sification schemes prone to the possibility of systematic
misclassification.
We attempt here to use time-series methods alone to
classify sources into general categories of burst-like (BL)
or stochastically variable (SV). These light-curve classes MD01
capturethevariabilitybehaviorofthetwomostcommon
MD02
extragalactic sources detected in image-differencing sur- g
veys, AGN and SNe. We present a novel method that MD03
r
separates BL and SV sources with high purity, using su- MD04
pervised machine-learning methods. Using multi-band s i
difference-flux in the gP1,rP1,iP1 and zP1 bands, we se- Field MD05 z
lgeacltaxByLhaonsdts.SVFofrrolmigh4t3cu61rvdeisffienreenacceh-imbaangde,swouerecsetsimwaitthe MD MD06 y
MD07
the fitnesses of analytical BL models relative to the OU
process, using both their estimated leave-out-one cross- MD08
validation likelihoods (LOOCV), and corrected-Akaike
MD09
informationcriteria(AICc). Weshowthattheuseofsim-
ple analytical models with suitably chosen priors, which MD10
mimic the approximate shapes of BL light curves (pre-
55000 55200 55400 55600 55800 56000
dominantly SNe), is sufficient for segregating them from MJD
SV, thereby obviating the need for exact models depen-
dent on specific BL subclasses. The model statistical Fig. 1.— The Pan-STARRS1 survey has a staggered 3-day
characterizations are combined across sources using a cadence in the gP1,rP1,iP1, and zP1 bands corresponding to 6
K-means clustering algorithm (Kanungo et al. 2002), to observations per month per filter, while yP1 is observed during
bright-time. The observations we use in this paper extend from
provide robust source classifications in each filter. The
2009September14till2011November17. InthispaperyP1isnot
filter-wise classifications are then averaged to give final usedduetotherelativelysparsecadenceascomparedtotheother
source classifications. Based on our BL and SV pho- filters.
tometric classifications and host galaxy offset cuts, we
define a photometric sample of SNe and AGN, which we The Pan-STARRS1 (PS1) telescope (Hodapp et al.
then use to define observational priors for future surveys 2004) is a 1.8 meter diameter telescope on the summit
such as LSST. of Haleakala, Hawaii with a f/4.4 primary mirror, and
This paper is organized as follows: In §2 we discuss a 0.9 m secondary, delivering an image with a diameter
the details of the survey and our data pipeline; in §3 we of 3.3 degrees onto 60, 4800×4800 pixel detectors, with
discuss time-series models used to describe BL and SV 10µm pixels that subtend 0.258” each (Tonry & Onaka
difference-flux light curves; in §4 we elaborate on com- 2009;Hodappetal.2004). Theobservationsareobtained
puting LOOCV and AICc to estimate model fitness; in through a set of 5 broadband filters g ,r ,i ,z ,y ,
P1 P1 P1 P1 P1
§5 we discuss our classification method, characterize the each with a limiting magnitude per nightly epoch of
3
∼23.5mag. AlthoughthefiltersystemforPS1hasmuch once the alert is registered, forced photometry can de-
in common with that used in previous surveys, such as tect both positive and negative fluxes. These criteria
the SDSS, there are substantial differences. For more remove the majority of “bogus” detections due to non-
technical details refer to Stubbs et al. (2010), and Tonry astrophysical sources, such as camera defects, cosmic
et al. (2012). rays, and image-differencing artifacts. The PS1 alerts
The PS1 survey has two operating modes, 1) the 3π are published to an online alerts database located in
survey which covers 3π square degrees at δ > −30 de- Harvard. Our automated pipeline then downloads the
grees in 5 bands with a cadence of 2 observations per alerts database to our local database servers at Uni-
filterina6monthperiod, and2)theMediumDeepSur- versity of Maryland on a nightly basis. The alerts are
vey (MDS) which obtains deeper multi-epoch images in thenprocessedandadditionalvalueaddedmeasurements
5 bands of 10 fields, each 8 square degrees, listed in Ta- are made on the data to enable easy characterization
ble1,designedforbothextensivetemporalcoverage,and of sources via a SQL-IDL-C++ pipeline (Fig. 2). The
full-survey stacked static-sky depth. Depending on the sources are automatically cross-matched with custom
weather, the accessible fields are observed with a stag- multi-band deep-stack catalogs (Heinis et al. in prep.)
gered 3-day cadence in each band during dark and gray toderivehostassociations. Otherstatisticssuchascolor
time (g ,r on the first day, i on the second day, evolution, and higher moments of magnitude and flux
P1 P1 P1
z onthethirdday,andthenrepeatwithg ,r ),and are also computed and stored in our database. Web
P1 P1 P1
in the y band during bright time. On average, the ca- pages that derive custom cuts on the data based on host
P1
dence(Fig.1)is6observationsperfilterpermonth,with properties, host offsets, color, magnitude, and time vari-
a 1 week gap during bright time, during which time the ability properties are also updated nightly. Our custom
Medium Deep fields are observed exclusively in y . We querypagecanbeusedtoquerythedatabaseanddisplay
P1
require the dense time-series (cadence≈few days) for ro- column-wise sortable results on a web page. The page
bustvariability-basedclassifications,therebymakingthe can also be used to visualize the data in our database
MD survey our survey of choice. using simple 2 dimensional plots or histograms that are
The PS1 MD data are processed using the image pro- created in IDL which are displayed on a web page. Fi-
cessing pipeline (IPP) located in Hawaii. The IPP per- nally, the transient alerts are classified based on their
forms flat-fielding and detrending on each of the indi- light curves using our time-series method discussed in
vidual images using white light flat-field images from a §5.
dome screen, in combination with an illumination cor-
rection obtained by rastering sources across the field of 2.1. Pre-Processing the Alerts for Classification
view. Bad pixel masks are applied, and carried forward Thetwomostcommonlyoccurringextragalactictime-
for use in the stacking stage. After determining an ini- varying sources are SNe and AGN, which fall under the
tial astrometric solution (Magnier et al. 2008), the flat- broad categories of burst-like and stochastically vary-
fielded images are then warped onto the tangent plane ing,respectively. However,thesebroadclassifications,in
of the sky, using a flux conserving algorithm. The im- combination with host galaxy offsets, also enable us to
age scale of the warped images is 0.250 arcsec/pixel. In discover more rare and exotic variables and transients.
the MD fields, all images from a given night are col- For example, nuclear BL sources should include tidal
lected with eight dithers. This allows the removal of de- disruption events (Gezari et al. 2012; Chornock et al.
fects like cosmic rays or satellite streaks, before they are 2014), off-nuclear BL sources may contain gamma-ray
combined into a nightly stack using a variance-weighted burst afterglows (Cenko et al. 2013; Singer et al. 2013),
scheme. Nightly stacks of images, each with a 8 square and off-nuclear SV sources may be offset AGN from a
degree field of view, as well as seasonal deep stack ref- post-merger recoiling supermassive black holes (Blecha
erence images are created, which are then transferred et al. 2011).
to the Harvard Faculty of Arts and Sciences Odyssey We identify extragalactic alerts by cross-matching the
Research Computing cluster, where they are processed 18,058 alerts detected in the first 2.5 years of the PS1
through a frame subtraction analysis using the photpipe MDS with galaxies detected in our custom multi-band
image differencing pipeline originally developed for the deep-stack star/galaxy catalogs (Heinis et al. in prep.).
SuperMACHO and ESSENCE surveys (Rest et al. 2005; Galaxies are detected in χ2 images (Szalay et al. 1999)
Gargetal.2007). Significantfluxexcursionsarethende- built from CFHT’s u-band and the five PS1 bands. The
tected in the difference images using a modified version detection threshold, defined by the χ2 distribution, is
of the DoPhot(Schechter et al. 1993; Rest et al. 2013) equivalent to a SNR of 1.9σ. The photometry is then
photometry package, and they are tagged as a source, if performed using SExtractor (Bertin & Arnouts 1996) in
they satisfy the following conditions: Kron elliptical apertures which are also used in cross
matchingalertswiththeobjectsinthecatalog. Thecat-
• Positive detections with a signal-to-noise ratio alogcontains≈107 objectswhichhavebeenclassifiedas
(SNR) ≥ 5 in at least three images within a time starsorgalaxieswithover90%accuracyforsourceswith
window of 15 days. magnitudes <24 mag, using an optimized SVM classifi-
cation scheme that takes into account the shape, color,
• Detections in at least two filters. andmagnitudesofthedetections(Heinisetal.inprep.).
Thusweonlyselectalertswithi <24mag,where
P1−host
• No previous alert at that position. the star/galaxy classification is reliable. We identified
4361extragalacticalertsusingthecatalog,whichwethen
While for source detection, DoPhot requires at least 3 characterized as SV or BL using multi-band difference-
consecutive positive detections within a 15 day period, fluxtime-series. Notethatwedonotincludeextragalac-
4
Fig. 2.—ThePS1-UMDdatapipeline. ThedataarerelayedfromtheIPPviathephotpipepipelinethatprovidestransientalerts,aswell
asperformsforcedphotometryandimagedifferencing. AtUMDthedataaredownloaded,enhancedwithstatisticalparameterizations,and
assimilated into SQL databases using a C++ framework, which are then queried using IDL or PHP for interactive web-based analysis.
tic alerts with unresolved hosts in our sample, such as the starting and ending points of light curves cannot be
quasars, or “hostless” alerts, those with either a faint subject to one of these criteria, we discard them after
host galaxy (i > 24 mag), or located outside the removing the erroneous difference flux points. Also, our
P1−host
elliptical region that defines their host galaxy. conditionwillnotweedouttheturn-onofBLlightcurves,
To characterize the sources, we model their difference- because only |y −y | < 10dy may not be satisfied
i i−1 i−1
fluxtime-seriesobtainedfromforcedphotometryineach due to a short rise-time, however, |y −y | < 10dy
i i+1 i+1
oftheg ,r ,i ,andz bands,andthencombinethe will be satisfied for all BL lightcurves, given the Pan-
P1 P1 P1 P1
characterizations across the filters. Forced photometry, STARRS1medium-deepsurveyobservationsrepeatonce
done by the photpipe pipeline, is obtained by performing every 3 days, in the g ,r ,i ,z bands. Finally, we
P1 P1 P1 P1
PSF fitting photometry on all difference images in each also transform the light curves such that the minimum
band,atthelocationofanytransientcandidate,inorder difference-flux is 0. This is done so as to make them
tofullyexploitallavailabledifference-fluxtime-seriesin- conform to the limits for the priors for BL light curves,
formationonthealert,usingdatafrompriortothealert especiallyforlightcurveswheretheBLsourcewasactive
detection. In our method, we decided to use difference- in the reference image.
fluxesinsteadofdifferentialmagnitudesbecausestochas- In our analysis we only use light curves which have at
ticlightcurvescanhavenegativedifference-fluxes(ifthey least n = 20 distinct difference-flux measurements post-
werebrighterinthereferenceimage),forwhichABmag- processinginanyfilter,soastoensurethatthisisatleast
nitudes cannot be defined. four times as large as the maximum number of parame-
Before we perform the classification on the difference- ters k used in any of the models (Maximum number
max
flux light curves, we pre-process them to remove arti- of parameters in any time-series model (§3) is 5). This
facts, as well as to make them conform with SV and BL is done to prevent over-fitting of the data, which may
model priors. Many difference-flux light curves also con- resultinmodelcomparisonsnotbeingmeaningful. Also,
tain singular large difference flux outliers which affect n = 20 is not a restrictive limit for classification pur-
model classification significantly. To remove them, each poses since this is a factor ≈2 smaller than the average
difference-fluxpointy inagivenlightcurveiscompared number of photometric measurements in any filter for
i
to its previous and next difference-flux values, y and all 4361 sources, which is ≈ 36. Fig. 3 shows the his-
i−1
y , and their photometric errors dy and dy ,with togram of number of distinct points in the g ,r ,i ,
i+1 i−1 i+1 P1 P1 P1
the criterion that at least one of |y −y | < 10dy , or z filters for all the PS1 transient alerts associated
i i−1 i−1 P1
or |y − y | < 10dy is satisfied for the difference- withgalaxyhoststhatpassourcuts. Inthenextsection
i i+1 i+1
flux point to be accepted as non erroneous. Since most we discuss the time-series models that we use to classify
difference-fluxerrorsaremuchlargerthanthiscutoffand the difference-flux light curves.
most differences in successive difference-flux values are
3. TIME-SERIESMODELS
much smaller than this, we can ensure that the light
curveisunaffected,whiletheoutliersareremoved. Since Since our goal is to classify extragalactic time-varying
sources into two broad classes, BL or SV, we assess the
5
MD01 MD02 MD03 MD04 MD05
1000 1000 1000 1000 1000
100 100 100 100 100
N N N N N
10 10 10 10 10
1 1 1 1 1
0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60
Number of Points Number of Points Number of Points Number of Points Number of Points
MD06 MD07 MD08 MD09
1000 1000 1000 1000
100 100 100 100
N N N N
10 10 10 10
1 1 1 1
0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60
Number of Points Number of Points Number of Points Number of Points
Fig. 3.— Sources plotted in the figure satisfy our criteria for classification: a light curve with n≥20 points in all 4 filters. MD10 has
nopointswhichsatisfyourcriteriaduetoonlyasinglefullseasonofcoverageinthefirsttwoandahalfyearsofthePS1MDS.
general shapes of the light curves by comparing their followed by a relatively slow decline in the flux in any
similarities to SN-like bursting behavior, or to AGN-like band (t ≈ 30 days). To resolve this, we employ a
fall
damped-random-walk type behavior. While fitting an Gamma distribution which is robust in modeling such
exactmodelinvolvesalargenumberofparameterswhich light curves (Fig. 4), reflecting varying degrees of asym-
maybeunknown,andmaynecessitatealargenumberof metrydependingontheshapek andscaleD parameters
data points, the general shape of a BL light curve can of the distribution. Another model is the Analytic-SN
be approximated to certain simpler analytical functional model, that uses distinct exponential rise and decline
forms (Gaussian, Gamma distribution, and generic ana- timescales, t and t , and is particularly well suited
rise fall
lytic SN model); and that of an SV light curve approx- to modeling non-Ia type SN light curves (Kessler et al.
imated by an OU process (Uhlenbeck & Ornstein 1930; 2010), although it is generic in its application. Despite
Andrae et al. 2013) as described in Table 1. In these thenon-specificnatureofthesemodels,wefindthattheir
models, we have ignored the effects of cosmological red- simple statistical descriptions of the difference-flux light
shift corrections and dust extinction. However this is curvesofBLsourcesaresufficientindistinguishingthem
acceptable since our goal here is to use the models only fromAGNwithlowcontamination. Theuseofthreedis-
todistinguishbetweencoherentsingle-bursttypebehav- tinct analytical BL models allows for a broader range of
ior from stochastic variability, while not assuming any BL light-curve shapes, and is comparable to using inde-
underlying physical processes for the sources. Since the pendentstatisticaldescriptionsofthelightcurvethrough
modelparametersspanseveralordersofmagnitude,their distinct parameterizations. Also, since the BL models
priors are chosen to be uniform in the logarithm. are compared with the SV model, only their relative fit-
The Gaussian is the simplest model that attempts to nesses in describing the data are important. Should the
model the overall flux from a BL source as the sum of necessity arise of classifying the objects into particular
a constant background α, and bursting behavior charac- sub-classes of the broader SV-BL distinction, or that of
terized by a Gaussian with amplitude β, center µ, and extracting particular details about the parameters of a
width σ. This however, does not account for the asym- source light curve, exact models for the sources (Kessler
metry in SN light curves; for example Type Ia SNe are et al. 2010; Kelly et al. 2010) must be included in the
betterapproximatedbyasharpriset ≈15days(Gai- comparisons, which although it is beyond the scope of
rise
tan2013;Ganeshalingametal.2011;Haydenetal.2010) this paper, is a direct natural extension.
6
0.5 100000
k=1 D=2 Gaussian -- OU Process --
k=2 D=2
0.4 k=3 D=2 80000 Gamma -- No Model --
k=5 D=1 s) Generic SN --
k=9 D=0.5 unt
o 60000
a (x) 0.3 ux (c
m Fl
m e 40000
Ga 0.2 nc
e
er
Diff 20000
0.1
0
0.0
0 5 10 15 20 0 100 200 300 400 500
x Time (days)
Fig. 4.—ArangeofBLlightcurveshapescanbemodeledusing Fig. 5.— An SN difference-flux light curve is reasonably well
Gamma distributions by varying the shape and scale parameters fitbyallthemodels,buttheBLmodelshavehigherLOOCVand
k,D. Thisisparticularlyapplicabletotheasymmetricriseandfall lower AICc as compared to the OU process resulting in the light
time-seriespatternsofSNlightcurves. curvebeingclassifiedasBL.
ThefluctuatingbehaviorofopticallightcurvesofAGN
is well described by an OU process (Kelly et al. 2009), a
first-order continuous-time auto-regressive process. The
processcanbedescribedintermsofadrivingnoisefield,
100000
parameterized by c, the square of the amplitude of the
Gaussian -- OU Process --
OU noise field, and a damping timescale τ (Bailer-Jones
2012). Mathematically, the evolution of the state vari- 80000 Gamma -- No Model --
able Z(t) of the OU process is given by the differential s) Generic SN --
nt
equation u
o 60000
c
x (
1 u
dZ(t)=c1/2dW(t)− (Z(t)−b)dt (1) Fl
τ ce 40000
n
whereW isaWienerprocess,andbisthemeanvalueof ere
theprocess. ForsimulatingtheOUprocessitself,weuse Diff 20000
the prescriptions from Bailer-Jones (2012). The method
uses posterior analysis to improve the estimate of the
state variable continuously, by using the observed flux 0
y at time-step t to compute the posterior distri-
k−1 k−1
0 100 200 300 400 500
bution of the state variable z , which is subsequently
k−1 Time (days)
used to update z . The OU process being a Gaussian-
k
Markovprocess, Z(t)ischaracterizedbyGaussianprob-
ability distribution function G(µ(Z(t)),V(Z(t))) where Fig. 6.—ExampleofanAGNlightcurvethatiswellfitbythe
µ,V are the mean and standard deviation of the state OUprocessandpoorlybytheBLmodels.
variable at time t. model likelihoods. Although for an OU process, the ac-
We also determine whether the light curves are well tual likelihood is computed differently from this (Bailer-
fittedbyaconstantmodelthatisrepresentativeofwhite Jones 2012), using only the photometric errors to com-
noise. In the event that none of the light curve models pute the likelihood for all models, is justified, since the
is significantly better than white noise, the light curves intent is to determine the model-mean that best mimics
are assumed to not pertain to any of the stochastic or the light curve shape. The probability P(yk|σk,θn) of
burstingcategories,andareclassifiedasNo-model(NM) observingadifferencefluxyk, assumingGaussianerrors,
sources. Figs. 5, 6, and 7 show examples of SN, AGN, is given by
andNMclassifiedsourcesandallthemodelfits. Ineach
case, the best models are chosen based on robust sta- (cid:18) 1 (cid:19) (f (θ )−y )2
tistical criteria, and the final source class decided using logP(yk|σk,θn)=log σ √2π − k n2σ2 k (2)
a clustering machine learning scheme described in the k k
following sections. wheref , y , andσ arethemodeldifference-flux, the
k k k
observed difference-flux, and the standard deviation es-
4. MODELLIKELIHOODANDFITNESSESTIMATION timates of the kth data point. For the OU process we
For all time-series models, including for the OU pro- use µ(Z(t )) or the mean light curve, in the place of f ,
k k
cess, we assume a Gaussian error model to compute the to evaluate its likelihood.
7
TABLE 1
Difference-flux models
Model Type Equation Parameters PriorDistributions
Gaussian BL Flux(t)=α+βe−(t−µ)2/σ2 α Ulog(0,105)(fluxcounts)
β Ulog(103,108)(fluxcounts)
µ Ulog(−102,104)(days)
σ Ulog(1,104)(days)
Gammadistribution BL Flux(t)=α+β(t−µ)k−1e−(t−µ)/D α Ulog(0,105)(fluxcounts)
DkΓ(k)
β Ulog(0,109)(fluxcounts)
µ Ulog(−102,104)(days)
D Ulog(1,102)
k Ulog(1,102)
Analytic-SNmodel BL Flux(t)=α+β e−(t−to)/tfall α Ulog(0,105)(fluxcounts)
1+e−(t−to)/trise
β Ulog(0,109)(fluxcounts)
to Ulog(−102,104)(days)
t Ulog(1,103)(days)
fall
trise Ulog(1,103)(days)
OUprocess SV dZ(t)=−1Z(t)dt+c1/2N(t;0,dt) τ Ulog(1,106)(days)
τ
whereZ hereisfluxcount. c Ulog(0,1014)(fluxcounts2)
b Ulog(0,108)(fluxcounts)
µ(Z) Ulog(0,108)(fluxcounts)
V(Z) Ulog(0,1014)(fluxcounts2)
No-Model WhiteNoise Flux(t)=C C Ulog(0,108)(fluxcounts)
20000 2k(k+1)
AICc=AIC+ (4)
Gaussian -- OU Process --
n−k−1
Gamma -- No Model --
15000
Generic SN -- The AICc is a correction to the AIC, that corrects for
s) the finite size of the dataset n relative to the number
nt
u 10000 of model parameters k. Note, that models that better
o
x (c represent the dataset have smaller AICc values. The
u LOOCV, another independent measure of model fitness,
Fl 5000
e is a measure of how well each difference-flux value can
c
n bepredictedusingtheremainingdifference-fluxdataand
e
er 0 hence, is a more robust statistical measure of model fit-
Diff ness as compared to the AICc. Re-stated, the cross-
-5000 validationcanalsobesaidtomeasurethepredictiveabil-
ityofamodelonagivendataset,andisanearlyunbiased
estimator(Kohavietal.1995)ofthemeansquarederror
-10000
0 100 200 300 400 500 for new observations. The LOOCV, more specifically, is
Time (days) theproductofthepiece-wiseprobabilityofobtainingin-
dividualdifference-fluxmeasurementsusingatime-series
model,whilesamplingtheparametersfromtheposterior
Fig. 7.—Exampleofadifference-fluxlightcurvethathasalarge constituted by the model over the remaining points in
numberofimage-differencingerrorsresultinginitbeingbestfitby
the time-series. In LOOCV estimation of a model over a
theNo-Model.
dataset y containing K points, the likelihood L of the
k k
kth data point is given by
To assess the fitness of the models, we estimate their
corrected Akaike information criteria (AICc) (Burnham
& Anderson 2002) and leave-out-one cross-validation 1 n(cid:88)=N
likelihoods (LOOCV) (Bailer-Jones 2012) over the Lk =P (yk|y−k,σ,η)≈ N P(yk|σk,θn,η) (5)
difference-flux data for each source, filter-wise. The AIC n=1
(Eq.3) is a quantification of the information lost when a
where η is the time-series model, σ is the error es-
model is used to represent a dataset. The AIC penalizes k
timate at each point, and θ are the model parameters
the maximum model log-likelihood lnL by a factor that n
drawninthenthiterationfromtheMarkovchainMonte-
depends on the number of model parameters k, thereby
Carlo sampling of the posterior probability distribution
accounting for over-parameterization of the dataset.
of the model over the other K −1 data points denoted
byy . TheLOOCVofthemodelcanthenbeobtained
−k
AIC =2k−2lnL (3) by multiplying the partition probabilities
8
Once the fitnesses of the models are estimated filter-
wise on the difference-fluxes using the AICc and the
k=K
(cid:89) LOOCV,weobtaintwoparameterspermodelfor5time-
LOOCV = L (6)
k series models in each of the four filters g ,r ,i , and
P1 P1 P1
k=1 z . First, we remove the NM best fit sources by com-
P1
Since the AICc and the LOOCV are measured inde- paring their model statistics with those of the BL and
pendentlyofeachother,theycanbeusedsimultaneously SV models. To do this we construct a relative sign vec-
to assess model likelihood, thereby reinforcing model fit- tor RV for each source, in each filter, using the AICc
i,f
ness assessment. The LOOCV for each model is esti- and the logarithms of the LOOCVs (which we designate
matedusingaMarkovchainMonte-Carlo(MCMC)using as LLOOCV):
a standard Metropolis-Hastings algorithm to sample the
posterior distributions (Hastings 1970). The model pa-
RV ={
rameters are sampled from known distributions and the i,f
posterior probability L p is evaluated, where p is the sgn(LLOOCV −LLOOCV ),
i i i Gaussian,i,f NM,i,f
prior probability and Li is the model likelihood in the sgn(LLOOCVGamma,i,f −LLOOCVNM,i,f),
ithiteration. Parametersforthei+1thiterationareac-
sgn(LLOOCV −LLOOCV ),
cepted with probability (L p )/(L p ), failing which Analytic,i,f NM,i,f
i+1 i+1 i i
the parameters from the ith iteration are retained. We sgn(AICcGaussian,i,f −AICcNM,i,f),
use a log-normal sampling distribution with a diagonal sgn(AICc −AICc ),
Gamma,i,f NM,i,f
covariance matrix, with σ2 = 10−4 uniformly across all
parameters, and all modelisi. We find that this choice of sgn(AICcAnalytic,i,f −AICcNM,i,f
sgn(LLOOCV −LLOOCV ),
a constant variance of the sampling distribution leads to OU,i,f NM,i,f
stable cross-validation likelihood values. In accordance sgn(AICc −AICc )
OU,i,f NM,i,f
with log-normal sampling requirements, the parameter } (7)
distributionsdefinedinTable1aretransformedbetween
−∞,∞ using a sigmoidal transform. where i is the object id, f is the filter, and sgn
The prior distributions for the BL and SV model pa- denotes the sign function, defined to be +1 for posi-
rameterscanbeassumedtobeuniforminthelogarithm, tive values and −1 for negative values. Ideally, for a
aswehaveinoursimulations,orcanbeobtainedbysam- BL source, RVBL = {+1,+1,+1,−1,−1,−1,±1,∓1}
plingtheparametersattheposteriormaximafortheBL since the BL models will have a larger LLOOCV, and
and SV models, for known SNe (BL) and AGN (SV) a smaller AICc when compared to the same for the NM,
training sets. The latter is advantageous if the entire while for an SV source the relative sign vector should
set of sources is well represented by the training set, in be RVSV = {±1,±1,±1,∓1,∓1,∓1,+1,−1}, i.e., the
that the number of iterations to convergence would be OU process better describes the light curve as com-
significantly reduced. However, we did not make this as- pared to any of the BL models or the SV models. For
sumption in order to allow for the classification of BL sources where the NM is the best model, RVNM =
and SV light curve types that may not occur in the veri- {−1,−1,−1,+1,+1,+1,−1,+1}.
fication set, and only took care to ensure that the limits We then compute RVi,f for all sources, which are ag-
ontheparameterrangessubsumedtheparametervalues gregated in filter-wise and fed into a K-means clustering
that could occur in the dataset. supervised-machine-learning algorithm, using the num-
Since the initial guesses for the model parameters in ber of centers K =3 in a swap method that is repeated
the MCMC may be far from the actual solution, a burn- over 100 iterations (Kanungo et al. 2002). The cluster-
in of 1000 iterations is employed for all model assess- ingalgorithmpartitionsthesourcesinthe8-dimensional
ments. We determined that a large number of burn-in RVi,f space, into Voronoi cells to determine the centers
iterationsisimportanttoensuresamplingnearthepeaks of the distributions for BL, SV, NM, by minimizing the
of the posterior distribution, and is particularly impor- sumofsquaresofthedistancesofpointsxl withincluster
tant while using prior parameter distributions that are Sm from the means of the clusters µm that correspond
uniform in the logarithm, as we have done here. We to the different classes of sources:
determined that 10000 post burn-in iterations were suf-
ficient for good model-fit convergence, after replicating
k
the results with 2000 burn-in iterations, and 20000 post (cid:88) (cid:88) ||x −µ ||2 (8)
l m
burn-in iterations. To ensure that the Markov chains
werewellmixed,wefurthercomputedtheGelman-Rubin l=1xl(cid:15)Sm
statistic (Gelman & Rubin 1992) over 10 parallel chains, EachsourceisthenassignedaclassC asBL,SV,or
i,f
foreachparameterineachmodel,andensuredthattheir NM depending on the center it is clustered around. The
valueswerelessthan1.1. ThecalculationoftheLOOCV squared-distanceofeachsourcepointiinfilterf,D =
i,f
is tedious and computationally expensive, and required |x −µ |2fromtheclusteringcenterµ isameasureof
i C,f C,f
us to parallelize our codes over a 300 core multi-node howreliablyitisclassifiedastheparticulartypeC,with
cluster. Inadditionourcodeswerewrittenground-upin adistanceofD =0beingthebest,andlargerdistances
i,f
C++ and optimized for quick run-time. The classifica- indicating less reliable classifications. C and D are
i,f i,f
tion of ≈ 7000 sources with ≈ 40 difference-flux points computed for each source, in each of the g ,r ,i ,
P1 P1 P1
in each of the four bands, required ≈4 hours. and z bands independently. We choose to classify the
P1
sourcesfilter-wise,andnotusingthestatisticalmeasures
5. CLASSIFICATIONMETHOD from all the filters at once, for the following reasons.
9
1. The behavior of each type of source, across the fil- modelcomparisonswhicharenotrelevanttoaparticular
ters, cannot be assumed to be uniform and hence, classificationtype,contributesignificantlytothenoisein
the clustering centers may differ significantly, theclusteringprocess. Wealsoattemptedaclusteringon
thedifferencesbetweenthemodelLLOOCVsandAICcs,
2. Clustering in some filters may be more noisy than
instead of on the signs of their differences, in a single
others, resulting in most sources being classified
clustering step, as well as in a hierarchical supervised
as no-model sources, thereby making these bands
method as discussed in this paper. However, we found
less favorable for classification purposes. For these
that in both cases, the number of misclassifications is
filters, the no-model center would be repeated in
larger due to the associated variance in the values of the
place of an SV or a BL center.
differencesinLLOOCVandAICc,whichismitigatedby
reducing them to binary statistics using the sign of their
3. Some filters may be less noisy and show clustering
differences alone.
only around two centers corresponding to BL and
Wecombinetheclusteringclassificationsfromthesec-
SV. Combining these filters with the noisier ones
ond clustering stage in each filter, by defining two mea-
resultsinboth,moreuncertaintyinclusteringclas-
sures; a quality factor C which is the average of classifi-
sification (larger D ) and a larger number of mis- i
i cations across the filters:
classifications. This is because, the uncertainty in
the clustering classification caused by one or more
(cid:80)
filters confounds the otherwise clear classifications C
from the others. As a result, the clustering cen- Ci = Nf i,f (10)
ters are poorly determined in the joint parameter filters
space of statistical parameters from all the filters. and, the average clustering square distance D across
i
By performing the clustering in each filter sepa- the filters:
rately, the classifications can be reinforced if they
show agreement across filters, and reflect the un- (cid:80)
D
vcearlutaeisn.ty otherwise, via smaller |Ci| and larger Di Di = Nffilteir,sf (11)
BLsourceshaveC closerto1whilestochasticallyvari-
Note, that it is favorable to assume more clustering i
able sources have C close to −1. D is a measure of the
centers in any filter than there are. For example, we i i
overallreliabilityoftheclassificationthatdecreaseswith
could assume that a certain filter has 3 clustering cen-
increasing D . Therefore, sources which are purely BL
ters corresponding to BL, SV, or NM while it may so i
will have C = 1,D = 0, while purely stochastic vari-
happen that one of the BL or SV centers is repeated, i i
ables will have C = −1,D = 0. Intermediate values
or two of the centers are relatively proximal, implying i i
of C indicate disagreements between some of the band-
that the clustering really occurs only around two cen- i
wise classifications, while larger values of D indicate a
ters. Post clustering, we filter out sources which have i
disagreement between the models in a given band.
been clustered around NM centers in at least 3 bands
they are detected in, and label them NM sources. The
5.1. Tests On a Verification Set
remaining sources then, are classified in at least 2 bands
they are detected in as either SV or BL. To detect their Totestourclassificationmethodweconstructedareli-
type more precisely than just using their comparisons to ableverificationsetwithadiverserangeofSNe(BL)and
the no-model, we construct another relative sign vector AGN (SV) in order to capture, as much as possible, the
BLSV comparing the fitness statistics of the BL mod- full range of their time-variability properties. For AGN,
i,f
els directly to those of the OU process for each source, we created a verification set from two sources: 1) 58
band-wise: UV-variability selected AGN with associations with PS1
alerts within 1(cid:48)(cid:48) from the GALEX Time Domain Survey
(TDS)(Gezari2013)withnoavailablespectroscopy,and
BLSV ={
i,f 2)125spectroscopicallyconfirmedAGNPS1alertsasso-
sgn(LLOOCVGaussian,i,f −LLOOCVOU,i,f), ciated with galaxy hosts from SDSS(Shen et al. 2011)
sgn(LLOOCV −LLOOCV ), and from a multipurpose Harvard/CfA program with
Gamma,i,f OU,i,f
the MMT to observe PS1 transients (PI Berger). The
sgn(LLOOCV −LLOOCV ),
Analytic,i,f OU,i,f GALEX AGN were selected from UV variability at the
sgn(AICc −AICc ),
Gaussian,i,f OU,i,f 5σ level in at least one epoch, and then classified us-
sgn(AICc −AICc ), ingacombinationofopticalhostcolorsandmorphology,
Gamma,i,f OU,i,f
sgn(AICc −AICc UV light curve characteristics, and matches to archival
Analytic,i,f OU,i,f
X-ray, and spectroscopic catalogs. The SN verification
} (9)
set consists of 131 spectroscopically confirmed Type-Ia,
We aggregate BLSV band-wise and perform a two- Type-Ib/c, Type-II, Type-IIn, andType-IIPSNefroma
i,f
center K-means clustering (K = 2) to segregate the BL combinationofPS1spectroscopicfollow-upprogramsus-
and SV sources. We find that this type of hierarchical ingGemini,Magellan,andMMTdescribedinRestetal.
supervised clustering, i.e. filtering out the NM sources (2013) and Berger et al. (in prep.). In order to test the
in the first stage and classifying the remaining sources performanceandefficiencyofouralgorithmintheclassi-
as BL or SV in the second stage, is more efficient as fication of AGN and SNe, we have constructed a diverse
compared to using all the model statistical comparisons and robust verification set that should be representative
concurrently in a single clustering step. This is because, of these populations in our sample.
10
(1+Density)% (1+Density)%
1 2 5 10 20 50 100 1 2 5 10 20 50 100
8 8
6 6
Di 4 Di 4
2 2
0 0
-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0
C C
i i
(a) (b)
Fig. 8.— Densities of verification set SNe (left) and AGN (right) on the Ci vs Di plane. Since AGN classifications for Di >4 occupy
boththeBL(Ci>0)andSV(Ci<0)regions,weonlyrelyonclassificationswithDi<=4. Asaresult,theSNeareclassifiedwith93.89%
completenessand90.97%puritywhiletheAGNareclassifiedwith57.92%completenessand95.00%purity.
Fig. 8 shows the contours of C vs D for the spectro- minimumoftheseasonalχ2sfortheverificationsetAGN
i i ν
scopicAGNandSNe. TheSNeclusteraroundtheregion andSNe. Ifχ2 =5.75isusedtoseparateAGNfromSNe,
ν
C ≥0.5 and D <4, while AGN predominantly occupy then 55.73% of the AGN can be recovered with 86.45%
i i
the regions defined by C ≤ 0 and D ≤ 8. In general, purity. We find that, using this value of χ2 to sepa-
i i ν
the degree to which an object is BL as opposed to SV rate AGN and SNe, results in maximal purity for the
increases with C . Some AGN light curves may show AGN sample. While this method yields a completeness
i
bursting-type behavior resulting in their being classified for AGN that is similar to that of our light-curve clas-
inmorethan1filterasBL,consequentlyhavingC >−1. sification algorithm, for SNe, the performance is much
i
Also, AGN clustering are less reliably classified as evi- worse, with 87.69% of SNe recovered at 47.22% purity.
denced by systematically larger D as compared to the This high contamination rate for SNe is due to the ex-
i
SNe. tensive overlap between the AGN and SN in the region
This is possibly because, the damped random walk χ2 <5.75. Hence,weconcludethatformaximumpurity,
ν
modelmaybeover-simplifyingthedescriptionoftheun- a more sophisticated method such as the one adopted in
derlying AGN variability. A more complex model such this paper is necessary.
as in Kelly et al. (2010) or Kelly et al. (2014), may be
required to capture the true variability. Another source 5.2. Final Classifications and Properties of
of confusion is that some of the AGN may resemble BL Extragalactic Sources
light curves if they fluctuate above the detection thresh-
We begin classifying our 4361 extragalactic transient
old for only a brief period during the observing baseline.
alerts by first selecting out sources which are clustered
Consequently, 47.88% of the verification set AGN clas-
around the NM center in at least 3 of the 4 bands. We
sified with D > 4, indicating unreliable classifications.
i find 570 such sources (NM sources hereafter). Visual in-
In order to maximize the purity of our classifications,
spection of the NM source light curves reveals that the
with a sacrifice to completeness, we use D = 4 as the
i majority are the result of noisy image-differencing light
bound for the classifications, below which 57.92% AGN
curves, most often due to large excursions in flux from
and 93.89% SNe, are recovered with 95.00% and 90.97%
imagedifferencingartifacts,andnotstatisticalerrorsdue
purities respectively. It is possible to include other pho-
tofaintfluxes, inmostofthebands. Thiscanbeseenin
tometric properties like color or host-galaxy properties
Fig. 9, which shows the minimum source magnitude in
to improve the completeness of the AGN classifications,
the g ,r ,i , and z bands for all sources, includ-
however, since the focus of our present work is to only P1 P1 P1 P1
ing that for NM sources (dark green), which barring the
usetime-variabilityasatoolforclassification,wereserve
brightest end of the magnitude distribution, follows the
this for future work.
overall magnitude distribution of extragalactic sources
In multi-epoch surveys such as Pan-STARRS1, it may
(black), indicating no strong biases toward fainter mag-
be possible to differentiate between AGN and SNe sim-
nitudes.
ply by comparing variability between observing seasons.
Fig. 13 shows C vs D contours for the 3791 extra-
For example, an SN will have only one season in which i i
galacticsourcesclassifiedSVandBL.Wedeterminethat
thereducedχ2 (χ2)islargerthan1, whileanAGNlight
ν there are 2262 SV sources and 1529 BL sources in the
curvecanshowvariabilityinbothseasonswithχ2 >>1.
ν dataset. Fig. 14 shows the distributions of SV, BL, and
In Fig.12 we test this simplified method by plotting the
NM by MD field, with SV being the most common class