Maximum-likelihood detection of sources among Poissonian noise

0.46 MB·English

by  I. M. Stewart

Astronomy&Astrophysicsmanuscriptno.1311papr January21,2009
(DOI:willbeinsertedbyhandlater)
Maximum-likelihood detection of sources among Poissonian
9
noise
0
0
2 I.M.Stewart1

n
a
J JodrellBankCentreforAstrophysics,UniversityofManchester,OxfordRoad,ManchesterM139PL,UnitedKingdom
e-mail:[email protected]
1
2
ReceivedJanuary0,0000;acceptedJanuary0,0000

]
M
Abstract. Amaximumlikelihood(ML)techniquefordetectingcompactsourcesinimagesofthex-rayskyisexamined.Such
I images,intherelativelylowexposureregimeaccessibletopresentx-rayobservatories,exhibitPoissoniannoiseatbackground
.
h fluxlevels.AvarietyofsourcedetectionmethodsarecomparedviaMonteCarlo,andtheMLdetectionmethodisshownto
p comparefavourablywiththeoptimized-linear-filter(OLF)methodwhenappliedtoasingleimage.Wheredetectionproceeds
- in parallel on several images made in different energy bands, the ML method is shown to have some practical advantages
o
whichmakeitsuperiortotheOLFmethod.SomecriticismsofMLarediscussed.Finally,apracticalmethodofestimatingthe
r
t sensitivityofMLdetectionispresented,andisshowntobealsoapplicabletosliding-boxsourcedetection.
s
a
[ Keywords.Methods:dataanalysis–Techniques:imageprocessing–X-rays:general

1
v
1. Introduction toussources(seeegDePonte&Primini1993;Dobrzyckietal
6
7 2000;alsotaskdocumentationfortheexsastaskdetectandthe
X-ray observatoriesfrom HEAO-2/Einstein onward have car-
2 SAStaskeboxdetect).
3 ried imaging cameras which have counted individual x-ray Thepresentpapershouldbeconsideredasanextensionto
. photons as they impinged upon a rectilinear grid of detector
1 paper I and also partly as an emendation. Several issues dis-
0 elements. The number of x-ray events detected per element cussed in paper I either had not at that time been completely
9 within a given time interval will be a random integer which
thoughtthrough,orwerenotdescribedasclearlyastheymight
0 followsaPoissonprobabilitydistribution.
havebeen.Theseissuesareitemizedasfollows:
:
v Asthesensitivityofdetectorsandtheeffectiveareasofx-
i ray telescopes have increased, so has the x-ray cosmic back- – Thedifferencebetweenthedistributionofsamplesofaran-
X
ground been increasingly resolved into compact sources. Eg dom function, and the distribution of its peak values, was
r
a a telescope such as XMM-Newton detects on average 70 notrealized.
serendipitoussourceswithin thefield ofview ofeach pointed – Theconceptofsensitivitywasnotadequatelydefined.
observation(Watsonetal2008).Afullcharacterisationofthe – The normalizationof the signal did notallow directcom-
x-ray backgroundwould shed light upon the evolution of the parisonbetweensingle-bandandmulti-banddetectionsen-
early universe(see Brandtand Hasinger 2005for a recentre- sitivities.
view of deep extragalactic x-ray surveys and their relation to – BothpaperIandthepresentpaperareconcernednotonly
cosmology). It is therefore of interest to find the best way to with numerical values of the sensitivity of various meth-
detect compact sources in these images - that is, the detec- ods of source detection, but also with practical, approxi-
tionmethodwhichoffersthebestdiscriminationbetweenreal matemethodsofcalculatingthesevalues‘inthefield’.One
sourcesandrandomfluctuationsofthebackground. may prefer a lengthy, complicated method for calculating
In a previous paper (Stewart 2006, hereinafter called pa- valuesfor a paper;such is unlikelyto be an equally satis-
per I) I discussed the general class of linear detection meth- factoryfieldmethod.InpaperItherewasinsufficientgrasp
ods, ie methods which form a weighted sum of the values of ofthisconceptualdistinction.
severaladjacentimagepixels,andoptionallyalsooverseveral – SensitivitiesfortheN-bandeboxdetectmethodwereincor-
different energy bands. A method for optimizing the weights rectlycalculated(seesection4.2foranexplanation).
was presented and shown to be substantially superior to the
sliding-boxmethodwhich,nodoubtbecauseofitssimplicity, The different signal normalization adopted in the present
hasbeenfrequentlyusedforcompilingcatalogsofserendipti- paper precludes direct comparison of the results herein with
those of paper I; however after appropriate rescaling the two
Sendoffprintrequeststo:I.M.Stewart sets of results are consistent. No numerical errors in paper I
2 I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise
werefound,despiteextensivereworkingofthesoftwareforthe theithvoxeldimensionstogivetherespectiveexpectationval-
presentpaper. uesB,S etc.Forthepurposesofthepresentpaperitismostly
i i,j
Theprimaryaimofthepresentpaperistopresentamethod assumedthatthecoordinatesrconsistofjustthespatial(x,y)
ofsourcedetectionbasedontheCashstatistic(Cash1979)and coordinatesofthe x-raydetector;in somesectionstheenergy
toshowthatthisisatleastassensitiveastheoptimizedlinear of the x-rayphotonsis added to the set. In this situation each
filter described in paper I. For this purpose, sensitivities of a voxelisboundedinthe(x,y)directionsofcoursebytheedges
variety of detection methods are compared across a range of of the physicalpixels of the detector. In the energydirection,
backgroundvaluesforbothsingle-bandandmulti-banddetec- voxelsareboundedbytheboundariesofanyenergyselection
tion scenarios. The paper also aims to demonstrate other ad- performedbytheuser.
vantages of the Cash source detection, and to discuss and in- The ultimate aim of source detection is to obtain the best
vestigatesomereporteddisadvantages. estimatespossibleofαand r forthemostnumerouspossible
0
Section2.1containsamathematicaldescriptionofthegen- subset of the M sources present. Satisfactory performance of
eralsourcedetectionproblem.Insection2.2a3-partdetection this dependsupon a couple of conditions. Firstly, the sources
algorithmisdescribed,inwhichamapofthedetectionproba- must not be confused - ie the density of ‘detectable’ sources
bilityateachimagepixelismadeasstepone,peaks(maxima) should be 1/∆r where ∆r is some characteristic size of
S S
≪
inthismaparelocatedinsteptwo,andasourceprofileisfitted S. For XMM-Newtonthisrequiresthe detectablesourceden-
toeachpeakinstepthree.Therelationbetweentherespective sitytobelessthanabout105deg 2.Thedeepestsurveystodate
−
probabilitydistributionsofthemapvspeakvaluesisalsodis- reachx-raysourcedensitiesofonlyatenthofthisvalue(Brandt
cussed.Detectionsensitivityisdefinedinsection2.3. andHasinger2005).Ontheotherhand,itwasshowninpaper
Cash-statisticsourcedetectionisdescribedinsection2.4.1. Ithatperturbationsofthebackgroundstatisticssetinatmuch
In section 2.4.2 it is shown how the Cash statistic may be shallowerlevelsofsensitivity;butstillalittlebeyondthedeep-
used at both stages 1 and 3 of the 3-part source detection al- estpracticalreachofXMM-Newton.Thereappearsthustobe
gorithm.Finally, in sections 2.4.3and 2.4.4, objectionsto the stillsomescopetoimprovesourcedetectioninthepresentob-
Cashmethodareaddressed. servatorieswithoutcomingup againstthe confusionlimit. To
Section 3 contains the results of the sensitivity compar- someextentalsoonecandisentangleconfusedsourcesduring
isons. In sections 3.1 and 3.2, details of the necessary calcu- thefittingprocess,bycomparinggoodness-of-fitforfittingone
lations are discussed. 1-band results are presented in section versustwo or more separate source profilesto a single detec-
3.3andN-bandresultsinsection3.4.Thedegradationofsome tion.
N-band sensitivities when the true source spectrum does not Ifwemayneglectconfusionthenwecandetectsourcesone
matchthatassumedisdiscussedinsection3.5. atatime.Themodelformulacanthusbesimplifiedto
An approximation useful for ‘in the field’ estimation of
e = B +αS(r r ). (1)
Cash sensitivities is presented in section 4.1. Finally, in sec- i i i− 0
tion 4.2, a deficiency of the sliding-box method is discussed, Thesecondconditionnecessarytosourcedetectionisthat
andamethodforremedyingitisproposed. theformofS is(i)known(atleastapproximately)and(ii)not
too similar to B. In x-ray detectors used to date, a significant
fraction of sources have angular sizes which are smaller than
2. Theory
theresolutionlimitofthetelescope.S forsuchsourcesisthen
2.1.Sourcedetection simplythePSFofthetelescope,whichcanbeestimatedapri-
ori,andwhichisoftenmuchmorecompact,iewithmorerapid
Inthesituationrelevantto thepresentpaperwe havea model
spatial variation, than the background.However, about8 per-
functione(r)definedoverthespacerby
centofsourcesdetectedbyXMM-Newtonforexamplearere-
M solvedorextended(Watsonetal2008).Detectionofaresolved
e(r)= B(r)+ α S (r r ) sourceismuchmoredifficult(oratleast,lessefficient)sinceits
j j 0,j
−
Xj S isdominatedbytheangulardistributionofsourceflux,which
isnotknownapriori.Inaddition,thereisofcoursenoupper
where B is called the background, each α is a scalar ampli-
limittothesizeofx-raysources,andinfactthedistinctionbe-
tudeand S isthe signalshape.We conceiveof a set of N po-
tweensourceandbackgroundislargelyaphilosophicalone.In
sitions r fori [1,N],eachwith an associatedmeasurement
i ∈ practice one has to impose an essentially arbitrary size cutoff
offluxc.Eachc isarandomintegerwithaPoissondistribu-
i i
and consider only variations in x-ray flux which have spatial
tion about a mean given by the value e of the model at that
i scales smaller than that cutoff to be sources, and take every-
position. From this comes the requirement that e 0 i. In
i
≥ ∀ thingelsetobebackground.
mathematicalshorthand
ModernCCDdetectorscanmeasurenotonlythepositions
M ofincidentx-rayphotonsbutalsotheirenergies.Thisopensan
c =e = B + α S (r r ).
h ii i i j j i− 0,j additionaldegreeoffreedomforS andBandthusthepossibil-
Xj ity of exploitingdifferencesin their respective spectra to sep-
In CCD x-ray detectors the c are obtained by accumulation aratethemwithevengreaterefficiency.Theaveragespectrum
i
withinarraysofvoxels,ievolumeswithinthephysicalcoordi- ofthesourcesinthe1XMMcatalogisquitedifferentfromthe
nates r. In thiscase the functions BandS are integratedover instrumentalbackgroundspectrumoftypicalinstruments(see
I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise 3
Carter and Read 2007 and Kuntz and Snowden 2008 for re- In the XMM-Newton data-product pipeline, both stages
cent data on XMM-Newton),and appearsalso to be different 1 and 2 are performed by the SAS task called eboxdetect,
to the spectrum of the cosmic background (see eg Gilli et al whereasthefinalstageisperformedbyemldetect.
2007). Howeverif itis true,as mostauthoritiesnowsuppose,
thatpracticallythewholeofthecosmicx-raybackgroundaten- Thefinaltaskmayperformsomeadditionalcalculations.It
ergieshigherthanafewkeViscomprisedofpointsources,one may calculate χ2 for the fit so as to aid in discriminatingreal
mightexpectthelatterdifferencetodiminishasfaintersources cosmicsourcesfromotherphenomena,suchas‘hot’CCDpix-
becomedetectable. els,whichmaygiverisetoexcesscounts(thisisnotatpresent
Parallel detection of sources in images in several energy donebyemldetect).Asisshowninsection2.4,useoftheCash
bandsofthesamepieceofskyiscomplicatedbythefactthat, statistic allows one to calculate a final value for the source-
althoughall(point)sourceshavethesamePSFatagivendetec- detectionlikelihoodafterthefit.
tor position and x-rayenergy,there is no correspondingcom-
monsource spectrum.Thismeansthatwe cannot‘locate’the This separation of detection and fitting stages is a useful
source S in the energydimension- we mustbe satisfied with procedure only where the image pixels are not significantly
locatingitspatiallyonthedetector. largerthanthe characteristicwidth∆r ofS. Thiscriterionis
S
Some of the detection algorithmsdescribedin the present infactachievablewiththeimagingcamerasofallthreecurrent
paperavoidmakinganyassumptionaboutthesourcespectrum; observatories primarily dedicated to detecting x-rays, namely
othersrequirethe user to choosea spectrum a priori. As will XMM-Newton, Suzaku and Chandra (Stru¨der et al 2001 and
beseenfromtheresults,thereisusuallyatrade-offinvolved,in Turner et al 2001 for XMM-Newton; Serlemitsos et al 2007
thatmethodswhichrequiretheuserto selectaspectrumtem- andKoyamaetal2007forSuzaku;andWeisskopfet al2002
plate may be better at detecting sources with spectra close to forChandra).
thattemplate,butworsewherethesourcespectrumisverydif-
ferentto the template;whereasa methodwhich makesno as-
sumptionswillusuallyperformwithaboutthesameefficiency
2.2.1. Peakdetection
onallsources.
Thepeak-detectionstagedoesnotseemtohavebeenthesub-
2.2.Likelihood-mapsourcedetection
jectofmuchstudy.Itisproblematicalfortworeasons.Firstly,
The present paper treats of a generic source-detection proce- the frequency distribution of values at the maxima of a ran-
durewhichfallsintothreeparts.Thefirstpartcalculatesamap dom function(call it the ‘peak’distribution) is not in general
givinganinitial,possiblycoarseestimateoftheprobability,for the same as the distribution of values over a set of randomly
eachpixel,thatasourceiscentredonthatpixel.Inotherwords, chosenpositions(whichwemightcallthe‘whole-map’distri-
ittests,foreach jthmappixel,thefollowingfluxmodel: bution).Ingeneralonewouldexpecttheformerdistributionto
besomewhatbroaderthanthelatter,sincetheprobabilitythat
e = B +αS(r r ). (2) a randompixelisa peakissmallforrelativelylowvaluesbut
i i i j
−
large for high values. Nakamoto et al (2000) give a method
ThistaskisaccomplishedbycalculatingsomestatisticU from tocalculatethedistributionofpeaksforarandomfunctionof
the count measurements ci taken from some subset of the N onedimension,butitisnotclearifthistheorycanbeextended
detectorpixels.Thespatialdistributionofthissubsetoughtnot to morethan one dimension.However,as is shownin section
tobemuchlargerthanthesizeofS -otherwiseoneisinclud- 2.4.4, the Cash maximum likelihood statistic does provide a
ingpixelswhichoneknowscontainnoinformationaboutS.In waytoestimatetheshapeofthepeakdistribution,althoughnot
other words, the pixels chosen to calculate the statistic at de- its normalization (which depends on the fraction of map pix-
tectorpixel jarebestlocatedwithinafieldofapproximatesize elswhicharepeaks).Thenormalizationmustbeestimatedby
∆rS,centredaboutpixel j. appropriatesimulations.
AftercalculatingthestatisticU forallpixels,thisfirstpart
of the detectionprocedureconvertsthis into a measure of the The second problem is that the definition of a peak in an
probabilitythatthevalueofU atpixel jwouldarisebychance arrayof samplesofa nominallycontinuousfunctionis some-
fluctuationofthebackgroundalone. whatarbitrary.Oneobviouschoiceofdefinitionistolabelany
Thesecondpartofthesource-detectionprocedureconsists jth pixelasa peakpixelif the valuec at thatpixelisgreater
j
ofapeakdetectionalgorithm.Theoutputfromthispartisalist thanthevaluesatallotherpixelswithinapatchofpixelssur-
ofpeaksintheprobabilitymap,whichareinterpretedassource roundingthe jth. Butwhatsize andshapeto makethepatch?
candidates. Itseemsclear thatthese oughtto matchthe size andshape of
The final task cycles over this list of candidates and at- S - too large a patch might miss additional, nearby sources;
temptstofitthesignalprofiletothefieldofc valuessurround- whereastoosmallrunsintotheoppositedangeroflistingmore
i
ingeach.Tobeofmaximumuse,thefittingprocedureshould thanonepeakforasinglesource.Also,itisnotclearhowthis
be capable both of separating a single candidate into two or choice mightaffect the frequencydistribution of peak values.
moreconfusedsources,andmergingmultiplecandidateswhich ThisdangerwaspointedoutbyMattoxetal(1996)inrespect
turnouttobelongtoasinglesource. ofsourcedetectioninEGRETdata.
4 I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise
2.3.Sensitivity ciated(throughthefunctionP(>U;0))withadefinitevalueof
U .
det
Theefficiencyofamethodofdetectingsourcesischaracterised
As an aside, note that in all cases in which the measured
byitssensitivity,which,broadlyspeaking,isthesmallestam-
data are randomintegers, p(U;α)is a sum of delta functions,
plitudeαwhichwillresultinasourcesignalαS beingdetected
thus P(>U;0) is a descending step function. Inversionof the
amongbackgroundof a given level. While this sentence con-
lattertoobtainU canclearlyonlygiveanapproximationin
det
veysthegeneralidea,aprecisemathematicaldefinitionisnec-
whichU issettothenexthighestvalueofU forwhich pis
det
essary before a comparison can be done between methods of
defined.Inthissituationonemayexpectthesensitivityofthe
sourcedetection.Asomewhatstiltedandobscuredefinitionof
methodto descendin stepswith decreasingbackground.This
sensitivity was presentedin paper I. The presentsection is an
effect can be seen in the plots of the sensitivity of the simple
efforttomakethisdefinitionclearerandmorerigorous.
sliding-boxdetectionschemeinfigures5,7and8.
Thesensitivityα ofadetectionmethodisheredefinedto
det
2.3.1. Definition bethatvalueofαwhich,forgivenBandS,givesrisetovalues
ofU suchthat U = U ,whereU isobtainedfromP as
det det det
Let us start by assuming that B and S are known a priori. described abovhe. Siee figure 1 for a diagrammatic example of
Formallyspeakingthisisnottrue:oneoughttomakeuseonly this.Thealgorithmforcalculatingsensitivityisthusasfollows:
oftheirbestestimatesBˆ andSˆ intheformulaforanysourcede-
tectionstatistic.Inpracticehowever,becauseBandS canoften 1. DecideavalueforP .
det
be quite well estimated by methods separate from the source 2. Invertthe relation P = P(>U;0)(solidlinein figure1)to
detection itself, it is convenient for the sake of simplicity to obtainU fromP ).
det det
maketheapproximationBˆ Betc. 3. An expressionforthe averagevalue ofU as a functionof
∼
EachofthedetectionmethodsdiscussedineitherpaperIor alphaisinvertedtogiveα =<U > 1 (U ).
det − det
inthepresentpapercanbeassociatedwithastatisticUwhichis
tobeevaluatedatpositionsofinterest:perhapsatthecentreof Note that other definitions are possible, for example one
eachimagepixel,orperhapsatasetofpositionsofcandidate could define α such that the median of U(α ) instead of
det det
sources. U at any position is evaluated, using a formula par- itsmeanwasadjustedsoastoequalU .Ihavepreferredthe
det
ticular to the detection method, from the measured counts c. definition in terms of the mean simply because it is generally
i
ClearlyU socalculatedisarandomvariablewhichwilloccur easiertocalculate.Howeverinthelimitofhighcounts,allthe
with some well-definedprobabilitydistribution p(U). We can varietiesofU herestudied tendto a gaussiandistribution,for
also define an integrated (more strictly, a reverse-cumulative) which the mean and median are equal. Within the range of
probabilityPsuchthat background values treated in the two papers, the median and
meandefinitionsempiricallyarefoundtoproducesimilarval-
P(>U)= ∞dU p(U). ues.
Z Themeandefinitionisequivalenttothe‘countsamplitude’
U
conceptemployedinpaperI;howeverthemediandefinitionis
InordertounderstandhowsensitivityrelatestoUitishelp- implicitinfigure12ofthatpaper.
ful to imagine an ensemble of measurementsof a source of a The curvesin figure 1 were constructed via Monte Carlo.
given amplitude α. Or one could simulate this situation via a At each iteration, simulated data values were generated for a
MonteCarloexperimentwhich,foreachmemberofitsensem- small (9 9) patch of pixels;the value of the appropriatede-
ble,generatesasetofcivalues.Theinputmodelforthisshould tection st×atistic U was then calculated for the patch. The end
consistofequation2andcanmakeuse oftheaprioriknown result of the Monte Carlo is an ensemble of U values. The
shapesofBandS.Foreachmemberoftheensemble,Uiscal- datavalueswerePoisson-distributedintegerswhosemeanval-
culatedfromthecivalues.TheresultingensembleofU values ueswere givenby a modelcomprisingflatbackgroundof1.0
will have a distribution which will in general be different for counts/pixelplusacentredsourceprofile.Thesolidcurveand
different a priori-assumed values of α. To reflect this, let us theerrorbarsrepresentMonteCarlo-generatedensembleswith
writeboththeprobabilitydensitypandthereverse-cumulative 106and105membersrespectively.
probability P(>U) with an additional α argument as p(U;α) Asmentionedinsection2.2,amulti-banddetectionmethod
andP(>U;α). whichdoesnotrequirethe userto makeanapriorichoiceof
We are interestedfirst ofallin the predicted P in the case source spectrum is expected to perform better at ‘across the
thatαequalszero-ie,whenthereisnosourceinthedetection board’detection,becauseofthefactthatcosmicx-raysources
field,onlybackground.Itisnecessarytoconsiderthisscenario don’t adhere to a common spectral template. Note however
becausetherandomnatureofthedatameansthat,atleastsome that,regardlessoftherequirementsofa(multi-band)detection
ofthetime,wewillunavoidablybelabellingasasourceacon- method, calculation of its sensitivity always requires the user
centration of counts which is merely a random fluctuation of to choose a spectral template. Eg in the Monte Carlo thought
thebackground. experiment described above, if α is non-zero, S needs to be
Thefirststepinsourcedetectionistodecidewhatfraction completelydefinedoverthewholespacer,whichinthemulti-
of these false positives we are prepared to put up with. This bandcaseincludesanenergydimension,inordertogeneratec
i
decisiondefinesacutoffprobabilityP whichinturnisasso- valuesandultimatelyavaluefor U .
det
h i
I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise 5
approximationforthe‘map’versionof thecumulativeproba-
bility function P exists for at least some of the statistics con-
sidered herein,whereas the ‘peak’ versionsof P must, for all
buttheCashstatistic,beestimatedviaalengthyMonteCarlo.
In addition, calculation of ‘peak’ distributions naturally also
requires the addition of a peak detector to the algorithm; the
MonteCarlomustcreatealargeimageratherthanasmallpatch
ofpixels;boundariesofsamehavetobeavoided,andsoforth.
Thethirdthingwhichmakesmapratherthanpeakcalcula-
tionsattractivearisesbecausethetruenumberdensityoffalse
detections(indeedofdetectionsofanykind)dependsnotjust
on the formula for U employed,but also on the nature of the
peak-detectionalgorithm.Havingtoconsidertheeffectofdif-
ferentalgorithmswithina‘space’withtwodegreesoffreedom
ratherthanoneisanunwelcomecomplication.
Noneofthesenicefeaturesofthe‘map’calculationwould
justify its use, however,if it were not for the fact, as pointed
Fig.1. The solid curve and the error bars (the latter shown in out by Nakamoto et al (2000), that the ‘map’ and ‘peak’ dis-
red in the electronic version) represent the results of separate tributionsofarandomfunctionbecomeasymptoticallysimilar
Monte Carlo experiments, as described in the text. The solid towards high values of the function. This convergence is ob-
curveshowsthereverse-cumulativeprobabilityP(>U;α)inthe served for example in figure 3. Since source detection by its
casethatthemodelconsistsentirelyofbackground(ie,α=0). naturedealswithextremevaluesofthenull-hypothesisproba-
The horizontaldotted line shows the chosen detection proba- bilityfunction,itwasfelttobe acceptabletoretainthe ‘map’
bilityP .Itsinterceptionwiththesolidcurvegivesthecutoff calculation for section 3, in which sensitivities of the various
det
valueU ofthedetectionstatisticU.Theerrorbarsshowthe methodsarecompared.Howeverwhenitcomestocalculation
det
approximate course of the probability density function p(U) of the absolute sensitivity of the Cash (or any other) statistic,
foracaseinwhichthemodelcontainsanon-zerosourcecom- ofcoursethetrue‘peak’versionofPmustbeemployed.
ponent.In the presentcase the modelsourceamplitude α has
beenadjusted such that U(α ) = U . This is the required
h det i det 2.4.TheCashstatistic
conditionforαtoequalthesensitivityα .
det
2.4.1. Definition
2.3.2. Mapvspeaksensitivity
Cash (1979) discussed a χ2-like statistic for Poisson data.
In paper I it was implicitly assumed that if one made a fre- Cash’s statistic was constructed in two stages. The first stage
quencyhistogram p (U) of U valuesat an ensemble of lo- comprisesthelog-likelihoodexpression
peak
cations of detected sources, that this would be related to the
N
probabilitydistributionp (U)ofU valuesatallimagepixels
map L= (c ln[e] e ln[c!])
bythefollowingproportionalrelation: i i − i− i
Xi=1
A
image
p (U)= p (U) where e is a modelfunction for the expectationvalue c . L
peak Abeam map can be sieen to consist of the logarithm of the producthofiithe
whereAimageistheimageareainpixels,andAbeamrepresentsa individualPoissonprobabilitydensitiesp(ci ei).Ifthemodele
|
constantknownasthe‘beamarea’,takentobeapproximately canbeexpressedintermsofoneormoreadjustableparameters,
equaltotheFWHMareainpixelsofthePSF.Ifthisweretrue, we can define the best fit values of these parametersas those
thecutoffprobabilityPdetcouldbeobtainedbymultiplyingthe which produce the largest value of L. Note that the factorial
maximumacceptablenumberoffalsedetectionsperunitimage termisnotafunctionofeandthuscanbeneglectedforfitting
areabythebeamarea.Inowrealizehoweverthatthisassump- purposes.
tionisnottrue:asdescribedinsection2.2.1anddemonstrated The second stage of Cash’s developmentwas to form the
infigure3,the shapesofthetwo distributionsaredifferent.It difference
followsthen thatone oughtto use the probabilitydistribution
ofU valuesatpeaksinthemap,ratherthanatallmappixels, U =2(max[L ] L ) (3)
Cash q null
−
incalculatingsensitivity.
This is a pity, because there are several attractive features where L is L calculated by setting all parameters of the
null
of calculating ‘map’ rather than ‘peak’ sensitivities. Three of model e to fixed, backgroundvalues and max(L ) for integer
q
these have present relevance. The first and most obvious of qisshorthandforthemaximumresultobtainablebyvaryinga
these‘attractions’istheretentionofconsistencywithpaperI. pre-selectednumberqofparametersofe.
Secondly, it is usually much easier and more straightfor- Solongasthenullhypothesismodelcanbeexpressedbya
wardtocalculate‘map’sensitivities.Anadequateclosed-form combinationoflegalvaluesoftheqfreeparameters,U 0.
Cash
≥
6 I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise
AccordingtoCash’stheory,underthisandotherconditionsdis- (xˆ ,yˆ ), for patch area A ; and let U (x,y) be the ideal,
0 0 peak ideal
cussed below, U in the asymptotic case is distributed ap- smoothfunction.Then
Cash
proximatelyasχ2,ieχ2forqdegreesoffreedom.
q
Just as with χ2 for gaussian data, the Cash statis- Uideal(xi,yi)= lim U(Amap)i,j
tic can be used not only to fit parameter values of the Amap→∞
model but also to obtain information about the probabil- and
ity distribution of those values. Cash for example uses it
to derive confidence intervals for fitted parameters. The U (xˆ ,yˆ )= lim U(xˆ ,yˆ ,A ).
ideal 0 0 0 0 peak
XMM-Newton SAS task emldetect (online description at Apeak→∞
http://xmm.vilspa.esa.es/sas/current/doc/emldetect.ps.gz) em-
Imagine the 3-part detection scenario of section 2.2 with
ploys the Cash statistic firstly to obtain best-fit values for the
theCashprescriptionusedbothtocalculatethemapofU val-
position, amplitude and optionally extent of sources in x-ray
ues in step 1 and to fit the source parameters in step 3. As
exposures,andsecondly,attheconclusionofthefittingroutine,
mentioned above, an advantage of the Cash statistic is that it
to calculate a detection probability for each so-fitted source.
can then be used to calculate a final value of U (thus of the
Calculation of the detection probabilityis done by setting the
null-hypothesis probability) for each fitted source. If the im-
sourceamplitudetozerointheexpressionforL .Theproba-
null age pixels are smaller than the characteristic size of the PSF
bilityofobtainingtheresultingvalueofU isthenequivalent
Cash (mentionedin section 2.2 as a conditionfor the usefulness of
totheprobabilityofobtainingthefittedvaluesofthesourcepa-
the3-partalgorithm),andifweassume(asseemsreasonable)
rameterswhenthereisnothinginthefieldbutbackground.
that the maxima in U are asymptotically paraboloid, then
ideal
DescriptionsoftheuseoftheCashstatistic forsourcede-
we expecttwo things:firstly, thatthe positionsoutputbystep
tection can be found in Boese and Doebereiner (2001) and
2,iepositionsofmaximainthearrayof‘map’valuesproduced
Mattoxetal(1996)togivejusttwoexamples.Theexpressions
bystep1,willbe,mostofthetime,within0.5pixelsofthefi-
developedinthepresentpaperarenotbasedoneithertreatment
nalfittedpositions;secondly,andbecauseofthis,thatthefinal
butcomparisonrevealsclosesimilarities.
U valueswillbeapproximatelythesameas,butslightlylarger
than,theU valuesatthestep2gridpositions.Theseresultsare
borneoutbyMonteCarloexperiments.
2.4.2. MapvspeakCash
Thesecondexpectation,plusthefactthat,asmentionedin
In the presentpaper,the Cash statistic is evaluated undertwo section2.3,therespectivetailsofthe‘map’and‘peak’distri-
differentscenarios:either at the centreof each image pixelto butionsofU areexpectedtobesimilar,togethermeanthatwe
makeamapofUCashsamples,oratthebest-fitlocationsofde- may investigate the sensitivity of Cash statistic source detec-
tectedsources.Inthefirstcase,thesourceistakentobecentred tion, orat leastcompareitwith othermethods,via the (much
onthepixelinquestion(equation2):theamplitudeistheonly moreconvenient)‘map’processalone.
fitted(iefree)parameter.Inthesecondcase,the xandycoor- Sincethis papershowsthattheCash statistic isatleast as
dinatesofthesourcecentrearefittedaswellastheamplitude. sensitive as, andmoreflexiblethan,anoptimizedlinearfilter,
It is as well to consider how the U values returned by these itmightbesupposedthata schemeofsourcedetectionwhich
respectivesituationsarerelated. employed Cash for stage 1 as well as stage 3 would be op-
Ineithercase,practicalconsiderationslimitthenumberof timum in practice. However, although it has been found con-
imagepixelsincludedinthecalculationofeachvalueofU venientto make use of this scheme as described abovein the
Cash
toarelativelysmallpatchsurroundingthepixelofinterest.As presentpaper,itcanbeshownthatthisisprobablynotthebest
can be seen from figure 6 however, the size of this patch has schemeforpracticalsourcedetection.Infactitisn’tnecessary
an asymptoticallydecreasingeffecton the resultof the calcu- tousetheCashstatisticforstep1aswellasstep3;noristhis
lation as this size becomes much larger than the characteris- likelytobeanoptimumuseofcomputerresources.Analterna-
tic size of the PSF. Consider then an ideal, continuous, two- tive methodis wellexemplifiedby the currentXMM-Newton
dimensionalfunctionU obtainedbycalculatingU ,with detection scheme, in which the detection statistic of step 1 is
ideal Cash
only the model amplitude free, at all points within the image calculatedusingasliding-boxformula(seeequationB.1).The
boundary.Samplesofthissmoothfunctiononthegridofpoints sliding-boxstatisticismuchquickertocalculatethantheCash
atthecentresoftheimagepixelsrepresenttheasymptoticval- statistic, and its relatively poor sensitivity is offset by setting
ues of the actual, practical map-U evaluation in the limit of the pass threshhold low enough to pass any possible source
large‘patch’size. The valuesatpeaksof U bearthe same through to step 3. Since step 3 does employ the Cash statis-
ideal
relationshipto the actuallycalculated valuesof peak-U atthe tic, the full sensitivity of the latter is retained without a great
endofthefinalfittingprocedure.Theadditionoftwomorefit- penaltyincomputingtime.
tedparametersdoesnotchangethevalueofU atthesepeaks: For purposes of the present paper, where the amplitude
onlytheappropriateprobabilitydistribution. aloneisfitted,thiswasdoneviaaNewton-Rhapsonmethod,as
Theseideascanbeexpressedmorecompactlyifweintro- detailed in appendix A. Where more than one parameter was
duce a little notation. Let U(A ) be the ‘map’ value of U tobefitted,theCashstatisticwasmaximizedviaaLevenberg-
map i,j
evaluatedatthepixelcentredon(x,y),forapatchareaA ; Marquardtmethod(followingthealgorithmdescribedbyPress
i i map
let U(xˆ ,yˆ ,A ) be the value at the fitted source location etal2003).
0 0 peak
I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise 7
2.4.3. ObjectionstoMLmeasuresofsignificance
Variousauthors(egFreemanetal1999, Protassovetal2002,
BergmannandRiisager2002,Pillaetal2005andLyons2007)
have criticized the use of Maximum Likelihood statistics for
thedetectionofsourcesorspectrallinesinPoisson-noisedata.
Thereseemtobetwomainconcerns:firstlythattheparameter
valueswhichreturnthenullhypothesislieontheboundaryof
thespaceofpermissablevalues;secondlythatthedatatowhich
aline(orsource)profileisfittedarealsothoseusedtoestimate
thebackground(nullhypothesis)signal.
The first objection appears to be a result of confusing
physicalwith statistical constraints. Statistically speaking,for
Poissondata,thereisnothingwrongwithafitwhichreturnsa
negativevalue for the amplitude,providedthe net modelflux
B + αS remains non-negative. Even if one specifies that the
physicalmodelatissuedoesnotallownegativeamplitudes(al-
thoughitiseasytoconceiveofsituationswhichcouldgiverise
toadipinfluxataparticularpositionorenergy),thefittedam- Fig.2. The error bars are the result of a Monte Carlo experi-
plitude αˆ and the model amplitude α are not the same thing: mentasdescribedinthetext.Theyshowthefrequencydistri-
theformerisonlyanestimatorforthelatter,thusisarandom butionofvaluesofthesourceamplitudeα,fittedbyuseofthe
variablewithacertainfrequencydistribution.Wherethemodel Cashstatistic.AllmembersoftheMonteCarloweregenerated
amplitudeissmallcomparedtothescatter inthefitted ampli- usingthesamemodelamplitude(3.0counts),indicatedinthis
tude,itisnecessarytoallownegativevaluesinorderforthefit- figure by the vertical dotted line. Negative values of fitted α
tedamplitudetobeanaccurateestimatorforthemodel.Figure were allowed; but the fit was not permitted to return a value
2demonstratesthis.(Figure2wasconstructedviaMonteCarlo less than the statistical limit (indicated with an arrow and the
inawaysimilartofigure1.Theensemblecomprised105mem- label‘Minimumα’),whichisdictatedbythenecessityforthe
bers.) total model flux to remain 0 for all pixels. The pile up of
≥
A negativevalueforαˆ whichis smallcomparedto itsun- theseboundedfitscanbeseeninthehighnumberswithinthe
certainty σ , where the physical model specifies that α 0, firstbin.Itcanbeclearlyseenthatitisnecessarytoretainthe
α
≥
merelyindicatesthatthemeasurementisconsistentwithzero. negativevaluesof fittedαif themeanofthe distributionisto
In this case one could use the sensitivity of the statistic as an beagoodapproximationtothemodelvalueofα.
upperlimit.Ifαˆ isnegativebutlargecomparedtoσ ,thecon-
α
clusion is rather that there is probablysomething wrong with
Thesituationisslightlymorecomplicatedinthemulti-band
the physical model (although, of course, such values will oc-
case. Because it is impossible to define a single source spec-
cur with a certain frequency by chance). In neither case are
trum,itmaybedesirabletoallowthefluxineachbandtobea
negativevaluesofαˆ somehowincorrect.Theonlyrequirement
freeparameter-ietofitthefluxα ineachkthbandseparately.
which the statistical theory imposes in cases where the noise k
Since one can easily conceive of a source which is bright in
is Poissonian is that the total flux (background plus source)
totalbutveryfaintinoneofmorespectralbands,itisentirely
maynotbe negative.Thusone canconclude,in contradiction
appropriatetoretainnegativevaluesforαˆ sinthefinalsource
toProtassovetal(2002),thatthenullvalueofα,iezero,does k
list.Asmentionedabove,allthatthisimpliesisthatthefluxof
notlieontheboundaryofthepermittedvaluesofαˆ.
thatsourceinthatbandisundetectable.
Preciselywhatonedoeswithnegativeαˆ valuesdependson
Clearly, a criterion to decide whether a given multi-band,
one’s aim. If that aim is to obtain the best estimate of α, as
free-spectrum detection should be retained or not should in-
in the above paragraphs,then negativeαˆs should be retained.
volve some functionof the fitted fluxes in all bands. It is just
Howeverif onedesiresto compilea listof candidatesources,
notquiteclearwhatthisfunctionshouldbe.Simplestmaybeto
thenfitsgivingnegativeαˆ oughtprobablytobediscardedfrom
simply addthe fluxesandretain onlythosesourcesforwhich
thelist.Innocaseisitcorrecttokeepthecandidatebuttamper
theresultispositive.
withtheαˆ value(by,forexample,settingittozero).
The validity of the second objection to maximum likeli-
Itshouldbepointedoutthatpolicytowardnegativeαˆ val- hoodfitting,namelythatthebackgroundandsignalamplitude
uesshoulddifferdependingonwhetherthedetectionisovera areobtainedfromthesamesetofdata,dependsonthecircum-
singlespectralbandoroverseveralinparallel.Forsingle-band stances of the detection. However, in XMM-Newton practice
detectionitisappropriatetodiscardcandidatesforwhichαˆ is at least, the fraction of image area occupied by identifiable
negative-indeed,ifadifferentalgorithmisusedforstage1of sourcesisusuallysmall.Inthiscaseitseemsareasonablepro-
the detection process, as in the XMM-Newton pipeline, such ceduretomakeanestimateofbackgroundfromthoseregions
candidates may never be submitted to Cash fitting in the first oftheimagewhichonehasconcludedprobablydon’tcontain
place. any detectable sources. The XMM-Newton pipeline performs
8 I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise
a simple 2-step iterative procedure to calculate a background
estimatebythismethod.Sinceoneisinthiscasenotusingthe
same data for both background estimation and source fitting,
thesecondobjectiondoesnotapply.
ProtassovetalproposethereplacementoftheMLstatistic
in these problemsby a Bayesian technique.While admittedly
more rigorous, their technique would require at least a fresh
Monte Carlo to be performedfor each source. Protassov et al
in their spectralline examplesdealin significanceson the or-
der of a few percent. To calibrate the distribution function of
their detection statistic, a Monte Carlo of a few hundred ele-
mentssuffices.Sourcedetectionhowevertypicallydemandsa
much more stringent significance level: 3 10 4 for 1XMM
−
×
(Watsonetal2003)forexample.MonteCarloensembleswith
on order 105 or more members are required to calibrate dis-
tributionsdowntothislevel.Sucha techniqueisnotpractical
forthecompilationofasourcecatalogwithupwardof50000
members.
Fig.3. Reverse-cumulativefrequencyhistograms from Monte
CarloensemblesoftheCashstatisticU .Continuouscurves
Cash
2.4.4. EmpiricalverificationofMLsourcedetection plotthetheoreticaldistributions,whicharejusttheχ2distribu-
n
tionsfortherespectivendegreesoffreedom,suitablynormal-
AnempiricallookatsomeMonteCarlodatashowsindeedthat
izedtotherespectivefractionofpixelsinvolved.Crosses,cir-
therawdistributionofvaluesoftheCashstatisticdoesnotap-
clesanddiamondsrepresentMonteCarlodata.Thesolidcurve
pear to conformvery closely to the ideal χ2 curve.This is so
(theory)andcrosses(MonteCarlo)showresultswhereonlythe
whetheronesetsthevalueofU tozerowheneverthefitted
Cash sourceamplitudeisfitted(=1degreeoffreedom),forasource
amplitude is negative (as in the examples of Protassov et al)
centredateachimagepixel.Thediamondsgivethedistribution
orwhetheronesimplyretainssuchvaluesunaltered.However,
of U valuesfor that small subsetof these ‘map’ pixelswhich
discardingnegative-amplitudefitsfromthepoolofsourcecan-
wereidentifiedaspeaks.Finally,thedottedcurve(theory)and
didates,asdescribedintheprecedingsection,bringsthedistri-
the circles (Monte Carlo) show the distribution of values ob-
butionoftheremainderveryclosetotheχ2 theoreticalcurve.
tainedbyfitting xandycoordinatesaswellasamplitude(=3
Figure3showstheresultsofsuchaMonteCarlo.
degreesoffreedom)foreachofthe‘map’peaks.
Forfigures3and4,itwasnotpossibletodealonlywitha
smallpatchofsimulated-datapixels,suchasforfigures1and
from a Monte Carlo. One can derive the number a posteriori
2, because in the presentcase it was necessary to performall
viaaplotsuchasfigure4.
threestepsofthegeneralizedsource-detectionalgorithm(sec-
Figure 4 was created as follows. An ensemble of (false)
tion2.2).Forthisreasonamuchlarger,circularfieldwasused,
sourcedetectionswas createdvia a Monte Carlo as described
witharadius(204pixels)similartothat( 210pixels)ofthe
∼ above in the description of figure 3. The detection software
field of view of the XMM-Newton EPIC cameras, when ex-
makes an estimate of the probability P that a source of that
pressedinthestandard4-arcsecpixelsofthedataproducts.No
amplitude or greater will occur by chance from background.
sourcesatallwereusedintheinputmodelfortheMonteCarlos
Ifthatestimateisproportionaltothetrueprobability,thenthe
of figures 3 and 4: all detectionsdepicted in these figuresare
proportionalitymaybeexpressedasfollows:
‘falsepositives’.
N(>P) M A f P
Thepenaltyforimprovingtheshapeofthedistributionby fields field peaks
∼
discarding negative fits is that one needs to estimate the ex- where N(> P) is the numberof sources in the ensemble hav-
pected fraction of such fits in order to normalizethe distribu- ing P or greater, M is the number of fields, A is the
fields field
tion. When considering the statistics of ‘map’-stage Cash (ie area of the field in pixels, and f is the fraction of image
peaks
fittingofamplitudealoneateachpixelofanimage),thisfrac- pixels which are identified as peaks (ie, positive-valued local
tionasymptotesto0.5inthelimitofhighbackground(seeeg extrema). (Note that the diamonds, circles and dotted line on
thezerointerceptofthesolidcurveinfigure3).Thisvaluehas figure3shouldalltrendto f asU approacheszero.)Taking
peaks
beenfoundto be an acceptableapproximationthroughoutthe logsandmanipulatinggives
rangeofbackgroundvaluesusedinthepresentpaper.
N L
Wherethesourcepositionaswellasamplitudeisfittedus- log + =log (f )
10 M A ! ln(10) 10 peaks
ing the Cash statistic, this simple approximationfor the frac- fields field
tionofimagepixelswhichyieldcandidatesnolongerapplies. wherethereplacementL= ln(P)hasbeenmadesincethisis
−
Inthiscaseonealsoneedstoknowthefractionofpixelswhich thevaluewrittenbyemldetecttothecolumnnamedDET ML
are peaks. There appears to be no way to estimate this apart intheoutputsourcelist.
I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise 9
‘falsepositives’.Thismayexplainthediscrepancywithpresent
results.
3. Sensitivitycomparison
3.1.Themodel
The flux modelused in the simulationwas based on equation
2.Thephysicalcoordinateshereofinterestarethetwospatial
coordinates x and y, representingpositionon the detector (or,
equivalently,position on a projectionplane tangentto the ce-
lestialsphere),plustheenergyEofthex-rayphotons.Samples
alongbothxandyaxesaretakentobeequallyspaced,located
indeed on an orthonormalarray, centred on the source centre
(x ,y ); but samples may be unevenly spaced in energy. It is
0 0
convenienttomodifyequation2tomakeuseofthreeindices,
onepercoordinate,viz:
e = B +αS .
i,j,k i,j,k i,j,k
Fig.4.Aplotwhichallowsthefractionofpeakpixels f to Bothiand jarerestrictedtotherange MtoMforsomesmall
peaks −
be estimated. This value is just the height to which the curve integerM.ThesymbolNishenceforthusedtodenotethenum-
asymptotesatsmallvaluesofL.Asdescribedinthetext,aflat berofenergybands.
slopeforsuchacurveindicatesthattheshape(butnotnecessar- ThesignalsamplesarerelatedtothePSFS(x,y,E)asfol-
ilythenormalization)oftheempiricalprobabilitydistribution lows:
agreeswith thatpredictedby theory.The blackbarsshowthe xi+∆x/2 yj+∆y/2
S = dxdyS(x,y,E )
r(eresdulitnsothfetheelebcetrsopnoikcevCerassihonfi)ttsihnogwsotfhtewraerseu;ltthseohfathlfe-tSonAeSbdaers- i,j,k Zxi−∆x/2 Zyj−∆y/2 k
where∆xand∆yaretherespectivesizesofthepixelsinthe x
tectionprocedure.Bothplotsaretruncatedsomewayshortof
andydirectionsandE =(E +E )/2,ietheaverageofthe
zerobecauseonlycandidatesourceshavinga‘map’probability k k,lo k,hi
boundingenergiesofbandk.S isnormalizedsuchthat
aboveagivencutoffweresubmittedtothefittingprocedure.A
numberof computationalreasonsled to the choice of a much 1 N
∞ ∞dxdy S(x,y,E )=1. (4)
largercutofffortheSASfittingtask(emldetect). ∆x∆y Z Z k
−∞ −∞ Xk=1
(Notethatthisisnotthesamenormalizationconditionasused
For each point on figure 4, N is calculated as the number for paper I.) Hence, α has units of x-ray counts and may be
ofsourceshavingthatvalueof Lorhigher.Theerrorbarsare interpreted directly as the expected total number of counts,
just √N.ThusN incrementsasoneproceedsfromthehighest summed over all N bands and over the entire detector plane,
Lvaluetothelowest. whichwillbegeneratedbythesourceS duringagivenexpo-
Resultsfortwosource-detectionproceduresareplottedon sureduration.
figure 4. The black bars show results from software written As for paper I, both ∆x and ∆y are set to the same value
for the present paper which performed all three stages of the of 4 arcsec, which is the same as the pixel size of the x-ray
source detection (section 2.2). The half-tone bars (red in the images producedas standard XMM-Newton products, and as
electronicversion)result fromuse of the XMM-NewtonSAS usedforsourcedetectioninboththe1XMMand2XMMcata-
taskseboxdetectfollowedbyemldetect.ActualXMM-Newton logs(Watsonetal2003and2008).
product file headers were used in the simulated image, expo- The PSF used, for both single-band and multi-band cal-
suremapandbackgroundmapfilesinordertofacilitateaccep- culations, was that correspondingto the pn camera of XMM-
tanceofthesebytheSASsourcedetectiontasks. Newtonontheopticalaxis,andataphotonenergyof1.25keV,
Theasymptoticflatnessofthefigure4plotsagainsupports whichisthecentreofband2ofthe1XMMcatalogx-raydata
theclosecorrespondencebetweentheactualdistributionofval- products(Watsonetal2003).InpaperIthe PSFwasallowed
uesoftheCashstatisticandtheχ2distributionpredictedbythe to vary in shape with energy band, as it is observed to do in
Cashtheory. XMM-Newton.Thisisdiscardedinthepresentpaperinorder
Finally, althoughit is beyondthe scope of the presentpa- toavoidconfusionbetweenthesmalleffectduetothechangein
pertoperformdetailedtestsonemldetect,onecanatleastsay shapeofS withenergyandthelargeeffectduetoitschanging
fromfigure4thatthevalueofDET MLisbeingcalculatedcor- magnitude.
rectly. Thisseems not to be consistent with the results shown The expected backgroundflux B in counts per pixel is
i,j,k
in figure 4 of Brunner et al 2008, which were also calculated taken throughout to be spatially invariant, but in general dif-
fromsimulateddata.However,theirsimulationwasmuchmore ferentwithineachenergyband.Sourceandbackgroundspec-
sophisticated:everyfieldcontainedverymanyfaintconfusing traused,plustheassociatedenergybanddefinitions,aretaken
sources, and a statistical procedurehad to be used to identify fromtable1ofpaperI.
10 I.M.Stewart:Maximum-likelihooddetectionofsourcesamongPoissoniannoise
3.2.Thesensitivitycalculation
Asnotedinsection2.3,thefirststep insourcedetection,thus
alsoincalculatingthesensitivityofthatdetection,istodecide
theacceptableprobabilityofafalsedetection.Here,asinpaper
I, a value of exp( 8.0)is used for this probabilitycutoff P .
det
−
This was the value used in the making of the 1XMM catalog
(Watsonetal2003).
ThenextstepistoinverttheappropriateexpressionforP(>
U;0)toobtainU .Inaminorityofcases,anaccurateclosed-
det
formexpressionforPisknown:inthesecases,Pwasinverted
byuseofRidders’method(Ridders1979)asdescribedinPress
etal(2003).Ridders’methodwaschosenratherthanaNewton-
typemethodbecauseofapaucityofexpressionsfor∂P/∂U.
In themajorityofcases, noclosed-formexpressionorap-
proximation for P is known which is accurate enough for
presentpurposes.Foreachdetectionmethodwithoutanaccu-
rateformulafor P, aMonteCarlotechniquewasusedto esti-
mateU .Anensembleofηsetsofc weregeneratedfroma
det i,j,k
modelcomprisingpurelythebackgroundBk.U wascalculated Fig.5. Minimum detectable source counts αdet, plotted as a
foreachsetofcountsvalues.ThefinalensembleofU values, function of background flux B; single-band case. The solid
m
for m [1,η], was sorted into increasing order. U is then line gives the results of the sliding-box technique; the dia-
det
∈
approximately that value of U for which (η m)/η P . monds show the results of detection by matched linear filter;
m det
− ∼
Clearlyηmustbechosensuchthatm 1. thecrossesshowtheresultsoftheCashmethod.
≫
Having obtainedU , the finalstep is to invertan expres-
det
sionfor U toderiveα .An accurateandinvertibleexpres-
det
h i
sion for U could be found for the majority of the detection there are many arrangementsof countsc which will produce
h i i
methods.Where not,a MonteCarlo techniquehadto beused the same U. The steps are still present for the other methods
for this step of the calculation as well. In this onerous and butareusuallytoofinelyspacedtobenoticeable.
somewhatproblematicprocedure,aRiddersiterationwasem- As one might expect, both the matched-filter and Cash
ployedtoconvergeonavalueofαdet.Foreachtrialvalueinthe methods are seen to be more sensitive (by a factor of at least
iteration sequence, an ensemble of sets of ci,j,k was obtained. 20%atanyvalueofbackground)thanthesliding-boxmethod.
Eachsetyieldedits valueofU. Themeanofthe ensembleof This is because they both make use of more information,
U values was tested against Udet and the trial α adjusted ac- namely the shape of S, to better discriminate between S and
cordinglyuntilconvergencewasjudgedtohavebeenreached.
B. What is not so expected is that the Cash method appears
tobemarginallybetterthanthematched-filtermethod.Atfirst
3.3.1-bandresults it seemedpossible thatthiswas an artifact,a resultofinsuffi-
cientconvergenceofthePowell’s-methodroutinewhichcalcu-
Threesingle-banddetectionmethodsarecomparedhere,viz: latestheweightsforthematched-filtermethod;however,tests
using much stricter convergence criteria failed to reveal any
1. Sliding-box:Uisobtainedsimplybysummingthec values
i significantdependence.Thematchedfilterdetectionisbydef-
withinthebox.
inition the optimum linear filter; however the Cash detection
2. Optimizedlinearfilter:Uisobtainedviaaweightedsumof
algorithmis based howeveron a nonlinear expression,which
thec,theweightsbeingchosentooptimizethesensitivity.
i isoptimizednotjustforthegeneralconditionsbutforthepar-
3. Cash-statistic detection: U given by the Cash statistic
ticulararrangementofcountswithineach(2M+1)by(2M+1)
whereamplitudeonlyisfitted.
detectionfield.Onemustassumethatthisaccountsforitsbetter
FormulasforP(>U;0)and U aregiveninappendixB.1. performance.
h i
Sensitivitiesofthethreemethodsarecomparedforarange Figure5iscomparableto figure6ofpaperI,with thead-
of values of background flux B in figure 5. All values were ditionofthepointsfortheCashresults.Somedifferencescan
calculatedforasquarearrayofpixelswithM =4. howeverbeobservedbetweenthesensitivity valuesdisplayed
The reason the sliding-box curve shows stepwise discon- infigure5andthosegivenintheearlierfigure.Therearethree
tinuities is because (as described in section 2.3) U for any reasonsfor this. Firstly, the normalizationexpressionfor S is
methodof detecting signals in Poissoniandata must be a dis- different:inpaperI,thenormalizationintegralforS extended
crete function, not defined for all values on the real line, due onlyasfarastheboundaryofthesquare(2M+1)by(2M+1)
to the discrete natureofthe countsvaluein anypixel.Allthe patchofpixelsinwhichcountswereconsidered;inthepresent
P functions, and any expression derived from them, thus de- papertheintegralextendsoverthewholeplane.
scendinstepfashion.Thisismostobviousforthesliding-box Secondly,asmoothapproximationtoP(>U ;0)wasused
box
method because U for this method is highly degenerate - ie inpaperI,ratherthanthetruestep-functionformusedhere.