Vi betraktar problemet med att uppskatta en rumsligt varierande densitetsfunktion, som motiveras av problem som uppstår vid storskalig radiologisk undersökning och anomalidetektering. I detta sammanhang är de densitetsfunktioner som ska uppskattas bakgrunden gamma-ray energispektra på platser spridda över ett stort geografiskt område, såsom kärnkraftsproduktion och avfallslagring platser, militära baser, medicinska anläggningar, universitet campus eller centrum av en stad. Flera utmaningar tillsammans gör detta till ett svårt problem. För det första kan spektraltätheten på en given rumslig plats ha både släta och icke-smooth-funktioner. För det andra är den rumsliga korrelationen i dessa densitetsfunktioner varken stationär eller lokal isotropisk. Slutligen, på vissa rumsliga platser, finns det mycket lite data. Vi presenterar en metod som kallas multiscale rumslig densitet utjämning som framgångsrikt tar itu med dessa utmaningar. Metoden är baserad på rekursiv dyadisk partition av provet utrymme, och därför delar mycket gemensamt med andra flerskaliga metoder, såsom wavelets och Pólya-träd tidigare. Vi beskriver en effektiv algoritm för att hitta en maximal a posteriori (MAP) uppskattning som utnyttjar senaste framsteg i konvex optimering för icke-smooth funktioner. Vi tillämpar multiscale rumslig densitet utjämning för verkliga data som samlats in på bakgrunden gamma-ray spektra på platser över ett stort universitet campus. Metoden uppvisar toppmodern prestanda för rumslig utjämning i densitetsuppskattning, och det leder till betydande förbättringar av effekt när den används i kombination med befintliga metoder för att upptäcka de typer av radiologiska avvikelser som kan få viktiga konsekvenser för folkhälsan och säkerheten. Nyckelord: radiologisk undersökning, täthetsuppskattning, rumslig statistik, Bayesian nonparametrics, totalvariation denoising, smält lasso Lost eller stulet radioaktiva källor utgör ett utmanande säkerhetsproblem. Brett används för industriell radiografi, sterilisering, och medicinsk avbildning, dessa källor är ofta dåligt säkrade (Gaffigan, 2012) och ibland stulna (Korshukin och Emery, 2006). För att förebygga farliga olyckor och upptäcka radiologiska spridningsanordningar (smutsiga bomber) innan de kan användas är säkerhetsbyråerna intresserade av att kontinuerligt övervaka stora områden för strålkällor. En enkel metod är att övervaka de totala strålningsnivåerna. Men dessa varierar naturligt i rymden, eftersom olika jord, sten och byggnadsmaterial kan innehålla mycket olika mängder naturligt förekommande radioaktiva material (NORM). Dessutom uppvisar olika detektorer mycket olika allmänna känslighet för strålning. Det är därför mer effektivt att övervaka energispektrumet för den detekterade strålningen istället, eftersom olika radioaktiva material avger gammastrålning vid distinkta energier. Men för att hitta en spektral anomali måste man först veta hur det normala bakgrundsspektrumet ser ut på alla rumsliga platser. Bakgrundsspektrumet är sannolikhetsfördelningen för energin från en gammastråle som avges av naturliga strålkällor på en given plats. Därför är ett viktigt statistiskt problem vid upptäckt av radiologiska anomalier behovet av att korrekt uppskatta en rumsligt varierande sannolikhetsfördelning för en uppsättning slumpmässiga variabler X (s), där s indexerar platsen. Vi presenterar en ny strategi för detta problem, kallas multiscale rumslig densitet utjämning, som utnyttjar en rekursiv dyadisk partition av provet utrymme X (s). Detta re-casts problemet med att uppskatta en densitet som problemet med att uppskatta en uppsättning binomial sannolikheter i samband med noderna i ett binärt träd. Vi inkluderar en prior på varje binomial sannolikhet som uppmuntrar rumslig utjämning, och vi beräknar en snabb MAP uppskattning under denna modell via en nyligen utvecklad algoritm för graf-baserad totalvariation denoising. Våra studier visar att metoden är snabb, mycket exakt för densiteter med svår flerskalig struktur, och skalbar till datamängder som är för stora för att lagras på en enda maskin. Vår huvudapplikation visar också att metoden leder till banbrytande prestanda när det gäller att upptäcka radiologiska anomalier, med hjälp av data om bakgrundsgammastrålning som samlats in under en sexveckorsperiod på ett stort universitetsområde. Våra data. Den dataanalys som beskrivs i detta dokument är en del av ett bredare arbete för att hjälpa dem som svarar först, som saknar förmåga att upptäcka och övervaka överdrivna strålningsnivåer, särskilt över stora geografiska områden. Nuvarande metoder är reaktiva, söker efter radioaktiva källor om de är kända för att gå förlorade, snarare än proaktiva, övervakning kontinuerligt för att upptäcka oväntade källor. I samarbete med lokala myndigheter utvecklar vi ett system för övervakning av strålning i komplexa stadsområden med hjälp av mobila 2 detektorer som används på befintlig infrastruktur, såsom sopbilar, stadsbussar, taxibilar eller postbärare. Målet är att utveckla en end-to-end system för övervakning och anomali upptäckt som inte kräver operatörsgränssnitt, en stor operatör utrustning börda, eller specialiserad utbildning. De data som analyseras i denna uppsats samlades in som en del av den första pilotstudien med detta system under juli och augusti 2012 vid University of Texas J.J. Pickle Research Campus (PRC), ett satellit campus som rymmer olika forskningsanläggningar vid University of Texas i Austin. Flera särdrag i Kina gör att bakgrundsspektrumet varierar geografiskt: flera tegelbyggnader har högre bakgrundshastigheter och det finns en radiologisk lagringsplats i nordvästra hörnet av campus. Genom att samla in data upprepade gånger under loppet av sex veckor, kan vi exakt kartlägga bakgrundsstrålning spektrum och observera dess förändring över tid, och upptäcka alla ovanliga spektral förändringar. För att samla in våra data använde vi en 2 × 2 tums cesiumjodid gammaspektrometer. Spektrometern spelade kontinuerligt in gammastrålarna, och svepte in dem i 4 096 diskreta energibehållare. 1 En laptop spelade in gammaspektrum, GPS-plats och tidsstämpel varannan sekund. Detektorn spelade in ungefär 30-50 gammastrålar per sekund. Detektorn kördes runt campus på en golfvagn en eller två gånger dagligen, med totalt 18 timmars observationer tagna över 41 enheter på 30 olika dagar. Se tillägg A för ytterligare uppgifter om data och förbehandling. Befintligt arbete med radiologisk undersökning. Förutom tillämpningen av lagstiftningen är radiologiska undersökningar viktiga i många andra säkerhets- och miljöscenarier. Inom medicinen kan de användas för att hitta en förlorad källa, t.ex. En medicinsk radioisotop som har stulits eller försvunnit från ett sjukhus. Vid gränssäkerhet kan de användas som en del av ett protokoll för att söka efter olagligt radiologiskt material på lastbilar vid en gränspassage eller containerfartyg i en hamn. Vid katastrofinsatser kan de användas för att bedöma strålningseffekterna av kärnkraftsolyckor som Tjernobyl eller Fukushima (Lyons och Colton, 2012). Tidigare källsökningsmetoder har fokuserat på att upptäcka radioaktiva källor med en mobil detektor när en tidigare radiologisk undersökning inte finns (Runkle et al., 2007; Jarman et al., 2010). Utan tidigare kunskap om bakgrundsstrålningens spektrum måste dessa system skilja mellan förändringar orsakade av mänskliga källor och naturliga variationer i bakgrunden, vilket minskar deras känslighet. Automatiserade system för övervakning av portaler, för användning vid gränsövergångar eller hamnar, har utvecklats (Runkle et al., 2009), och portalövervakare är nu allmänt utplacerade. Men även om vissa städer har genomfört helikopterbaserade strålningsundersökningar för katastrofberedskap (Wasiolek, 2007), dessa är opraktiskt dyra för långsiktig strålningsövervakning. Tidigare arbete på långsiktig kartläggning och anomali upptäckt inkluderar Aucott et al. (2013), som försökte kartlägga bakgrunden och korrelera spektral funktioner med 1 Eftersom energibehållare är diskreta, vi inte uppskattar en densitet i strikt mening, bara en mycket fin histogram. För att undvika onödig komplexitet hänvisar vi dock till "täthet" och "densitetsutjämning". 3 funktioner som upptäcks av kameror och LIDAR, såsom närliggande byggnader och byggmaterial; Reinhart et al. (2014), som använde spektraljämförelseförhållanden för att upptäcka spektralförändringar i stora 250 × 250 meter rumsliga celler, och Reinhart et al. (2015), som använde Kolmogorov-Smirnov tester för att uppnå högre effekt med mycket mindre rumsliga celler, vilket ger en bättre förmåga att upptäcka och lokalisera källor. De senare två metoderna lånar inte information över rymden: varje cell behandlas som oberoende av de andra, så varje cell kräver tillräckliga bakgrundsobservationer för att få en bra bakgrundsspektrumbedömning. Syftet med denna artikel är att visa hur en ny metod kan ta itu med denna brist och förbättra övergripande detektionsprestanda. Flera faktorer gör vårt problem svårt. För det första, bakgrundsspektra vi försöker uppskatta har både släta och spetsiga funktioner. De beskrivs inte väl av någon enkel klass av densitetsfunktioner eller av en enda bandbreddsparameter. Utanför ett rumsligt sammanhang behandlas detta ofta genom lokalt adaptiv densitetsuppskattning (t.ex. Sain och Scott, 1996). Det är ett svårt problem i sig, oberoende av den komplikation som den rumsliga variationen medför. För det andra är den rumsliga variationen i bakgrundsspektra icke-stationär och lokal anisotropisk. En anledning till detta är att olika byggmaterial har olika kombinationer av naturligt radioaktiva isotoper, och det finns stor variation i byggmaterial över platser: vissa radioaktiva isotoper är vanliga, medan andra är sällsynta (Ryan et al., 2014). Slutligen finns det på vissa platser mycket få uppgifter, eftersom systemet för insamling av bakgrundsuppgifter inte ofta besökte dessa regioner. Denna databrist innebär att en noggrann rumslig interpolering av bakgrundsspektrumet ofta är nödvändig för att upptäcka avvikelser på sällan besökta platser. Vår föreslagna flerskaliga metod för densitetsutjämning tar framgångsrikt itu med dessa utmaningar. Konceptuellt innebär det tre steg: 1. Dela upp det övergripande problemet i delproblem genom att rekursivt dela upp provutrymmet i en serie halvrymder, var och en beskriven med en villkorlig sannolikhet. 2. Utgångspunkten är följande: Jämna halvrymdsannolikheter genom att beräkna en snabb MAP uppskattning under en tidigare som uppmuntrar rumsliga grannar att vara liknande. 3. Vad är det som händer? Sammanfoga de utjämnade halvrymdssannolikheterna för att ge skattningarf (s) vid varje plats. Denna "splitt/smooth/smooth" struktur minskar täthetsutjämningsproblemet till en uppsättning av separata skalor utjämningsproblem, som var och en kan lösas genom konvex optimering. Den är därför extremt snabb och skalbar till mycket stora datamängder, och kan utnyttja två kompletterande strategier för parallellisering. För det första kan data för olika rumsliga platser lagras på olika maskiner, med minimal kommunikation mellan dem. För det andra kan olika scallar utjämningsproblem, motsvarande olika halvrymdsannolikheter, lösas på olika maskiner och kombineras enkelt för att ge ett slutligt svar. Vi figur 1 : Två exempel rumsliga platser (cell 14, översta raden; och cell 62, nedre raden) för Kina data som visar jämförelsen mellan den globala densitet uppskattningen och den lokala densiteten uppskattningar med hjälp av vårt förfarande för multiscale rumslig densitet utjämning. De två vänstra panelerna jämför träningsdata i varje cell (grå histogram), den lokala densitetsuppskattningen baserad på flerskalig rumslig utjämning (solid grå kurva) och den rumsligt invarianta globala densitetsberäkningen (dalerad svart kurva). (Dessa densiteter omskalas för att visa den förväntade talhastigheten.) De två mellersta panelerna visar skillnaden mellan de globala och lokala uppskattningarna nära den breda toppen i spektra mellan kanalerna 50 och 200. De högra två panelerna visar platserna för varje plats på en större karta över campus. 5 utnyttja båda källorna till parallellism i vår analys. Med dessa besparingar blir det möjligt att beräkna enorma radiologiska datamängder. Figur 1 ger en ögonblicksbild av metodens prestanda på våra data, samt viss intuition om de statistiska utmaningar vi har nämnt. Figuren visar täthetsrekonstruktionen för två platser (märkta 14 och 62) för vilka en relativt stor mängd data finns tillgänglig (eftersom kartan har delats upp i ett mycket grovt rutnät). De råa gammastrålarna för dessa två platser uppvisar synbart olika bakgrunder nära den breda topp i spektrumet som uppstår mellan energikanalerna 50 och 200. Vår procedur rekonstruerar troget dessa skillnader utan den alltför stora variationen i de råa histogrammen, men ändå lyckas den fånga andra "spiky" drag i densiteten utan översmoothing. Anslutning till befintligt arbete. Vår metod bygger på två stora, oberoende organ av statistisk litteratur: densitetsuppskattning och rumslig utjämning. Det finns en betydande mängd arbete på villkorlig densitetsuppskattning, av vilken rumsligt varierande densitetsuppskattning är en delmängd. Det vanligaste arbetet med detta ämne är baserat på kärnmetoder (se t.ex. Bashtannyk och Hyndman, 2001, och hänvisningarna där). Dessa tekniker innebär dock vanligtvis att uppskatta en densitet som förändras smidigt som en funktion av kovariater. De är inte direkt anpassade till det rumsliga sammanhanget, och i vilket fall som helst är den rumsliga variationen i våra data inte jämn. Dessutom gör traditionella (icke-adaptiva) kärnmetoder dåligt vid uppskattning av densiteter med flerskalig struktur som de vi möter, som innehåller både spiskiga och smidiga funktioner, både i rumsliga och energidomäner. Ingen felfri kärnmetod verkar lämplig för vårt problem. Det finns också en betydande kropp av Bayesian arbete på villkorlig densitet uppskattning med hjälp av den beroende Dirichlet processen och dess variationer (MacEachern, 2000). Men dessa metoder expanderar dåligt, vilket gör dem oinskränkt långsamma för vårt problem. Tekniken att representera en sannolikhetsfördelning via en rekursiv dyadisk partition har använts i stor utsträckning, till exempel i flerskaliga modeller för Poisson och multinomiell uppskattning (Fryzlewicz och Nason, 2004; Jansen, 2006; Willett och Nowak, 2007), och i icke-parametrisk Bayesian slutsats via Pólya-träd tidigare (t.ex. Mauldin et al., 1992; Lavine, 1992 Lavine, 1994. Konceptuellt är vår uppsats närmast relaterad till förslaget av Zhou och Hanson (2011), som par ett Polya-träd före med en villkorlig autoregressiv (CAR) innan att uppskatta en rumsligt varierande överlevnadsmodell. Men deras är en helt Bayesian modell och måste beräknas av Markov-kedjan Monte Carlo. Detta skalar dåligt till mycket stora rumsliga gitter av det slag som krävs för att karakterisera finskaliga radiologiska funktioner över en hel stad. Vår metod är inte helt Bayesian, och den passar genom konvex optimering, som skalar mycket mer fördelaktigt. Dessutom, på experiment som är tillräckligt små för att metoden Zhou och Hanson (2011) ska passa, vår metod har jämförbar eller bättre noggrannhet, men är storleksordningar snabbare. Det pris vi betalar för skalbarhet är att, till skillnad från en helt Bayes strategi, vår metod inte ger fel bars-bara en punkt uppskattning. 6 Spatial utjämning är ett stort och moget fält, och vi försöker inte en omfattande översyn. Vårt tillvägagångssätt är nära relaterat till en teknik som kallas totalvariation denoising från bildbehandling litteratur (Rudin et al., 1992), som också är känd i statistik som smält lasso (Tibshirani et al., 2005). Vi drar stor nytta av det senaste arbetet om beräkningseffektivitet för denna klass av optimeringsproblem, inklusive Tibshirani och Taylor (2011) (2014), och Tansey och Scott (2015). Speciellt använder vi algoritmen från Tansey och Scott (2015) (som i sig är starkt motiverad av diskussionen i Wang et al., 2014) för att lösa en rad optimeringsproblem som kombinerar en binomial sannolikhet tillsammans med en totalvariation straff över noderna i en oriktad graf. Vårt tillvägagångssätt är också begreppsmässigt relaterat till Bayesian inference för areal data modeller. Specifikt liknar det villkorliga auto-regressiva (CAR) modeller (Besag, 1974), som också påverkar rumslig utjämning genom att avskräcka stora parvisa skillnader över kanter i en graf. Icke-parametriska Bayesianska metoder för rumslig utjämning undersöktes av Gelfand et al. (2005), Reich och Fuentes (2007) och Rodriguez m.fl. (2010), bland många andra, och vi hänvisar intresserade läsare till dessa dokument, och deras referenser, för mer detaljer. Avsnitt 2 beskriver den övergripande strukturen för vår "split/smooth/smooth" strategi. Avsnitt 3 beskriver vår strategi för rumslig utjämning och vår passande algoritm. I avsnitt 4 presenteras simuleringar som visar att vår metod är mycket effektiv när det gäller att uppskatta en rumsligt varierande densitetsfunktion. Avsnitt 5 beskriver vår huvudsakliga tillämpning vid upptäckt av radiologiska anomalier. Avsnitt 6 avslutas med några avslutande kommentarer. Översikt av förhållningssätt I detta papper antas utrymmet vara diskret och representeras av en oriktad graf G = (V, E ), där en kant (r, s) innebär att platser r och s är grannar. Vi använder p =  och V för att beteckna antalet hörn och d =  och antalet kanter. Låt f (s) vara uppsättningen av densiteter över platser, alla med gemensamt urval utrymme B. På varje plats observerar vi data En nyckelenhet i vårt papper är användningen av en rekursiv dyadisk partition av provet utrymme B. Vi rekursivt konstruera nivå-k partition av B, betecknad Π (k), via en bijection med alla längd-k binära sekvenser γ  och {0, 1} k, enligt följande. Låt partitionen nivå-1 som Med tanke på partitionen på nivå k, definieras partitionen på nivå k + 1 genom att för varje γ  på 0, 1} k ange paret (B γ0, B γ1 ) som uppfyller B γ0 B γ1 = B γ och B γ0 B γ1 = ∅. Här γ0 (γ1) är ny binär sekvens definieras genom att lägga till en 0 (1) till slutet av γ. Vi hänvisar till Bγ (eller bara γ) som en förälder nod, B γ0 som vänster barn, och B γ1 som rätt barn. För att ge ett konkret exempel, om B är enhetsintervallet (dvs. Nivå-0 partitionen), nivå-1 partitionen kan vara Π ( Tillsammans med en uppdelning regel för att bygga barn noder B γ0 och B γ1 från B γ, definierar denna rekursiva process ett djup-K binärt träd vars kompletta uppsättning noder indexeras av alla binära sekvenser γ längd inte mer än K. Det finns många möjliga delningsregler; vi väljer B γ0 och B γ1 att ha lika längd, förutsatt att B är kompakt. Detta träd B (K) möjliggör en "split/smooth/smooth" uppskattningsstrategi, enligt följande. För närvarande tappar vi det rumsliga indexet s, och betraktar densiteten f på en enda plats (t.ex. Panel A i figur 2. Antag att X  till f, och att B γ är några (icke-terminal) nod i B (K). Vi fokuserar på den villkorliga sannolikheten att X faller i vänster-barn nod B γ0, med tanke på att det är i B γ : Vi uttrycker sedan sannolikheterna vid slutnoder av trädet som produkten av w γ s som man passerar trädet börjar vid roten. Dessa terminal-node sannolikheter bildar en mycket fin diskret approximation till f. Nu anta att x = (x 1,. .., x n ) är en vektor för n i. i.d. Prover från f, och låt y (x; b) = #{x i b} vara antalet av dessa prover som faller i uppsättningen b  på B. Om B γ är en icke-terminal nod, då tydligt där y γ = y(x, B γ ) för kort. På så sätt har vi reducerat problemet med densitetsskattning till en uppsättning oberoende binomiala problem, en för varje nod i B (K). Avgörande är att varje w γ kan uppskattas oberoende av varandra. Detta förklarar de "splitta" och "sammanfoga" aspekterna av vår uppskattningsstrategi. I delsteget bildar vi separata binomial-skattningsproblem för varje w γ. I sammanfogningen steg, efter att ha uppskattat w γ s, vi rekursivt beräkna terminal-node sannolikheter från den villkorliga uppdelning sannolikheter. Den saknade ingrediensen här är rumslig utjämning. För att förklara detta, vi återinföra det rumsliga indexet s och överväga insamlingen av dela sannolikheter w γ för alla träd noder γ och platser s  på V. Panel B i figur 2 visar en leksak exempel för noden γ = 1, där den underliggande densiteten skiljer sig för varje plats i en 2 × 2 rutnät. Konkret, på plats 3, har vi y 10 = 421 och y 11 = 127. Vår metod tillämpar rumslig utjämning till batchen av dela sannolikheter W γ = w (s) γ : s  och V associerade med en enda splitting nod γ, över alla rumsliga platser, före "sammanfogning" steget. Vi gör detta med hjälp av den binomial sannolikhet definieras av Equation (1) tillsammans med en strafffunktion som uppmuntrar rumslig jämnhet. Resten av dokumentet kommer att beskriva detta steg i detalj. 3 Spatial utjämning via graf-baserad denoisering Genom konstruktion, batchen av dela sannolikheter för en trädknuta, w (s) γ : s  på V, kan uppskattas separat från batchen för en annan trädknuta. Därför anser vi att utjämningsproblemet för en enda icke-terminal nod B γ  på B, och vi tappar indexet γ. Låt y (s) y (x (s), B γ0 ) vara antalet i vänster barn på plats s, och låt m (s) y (x (s), B γ ) vara antalet i föräldern. Antalet vänsterbarn följer en binomial fördelning, med tanke på föräldratalet: Målet är att uppskatta den villkorliga vänstersplit sannolikhet w (s) på varje plats. Vi tar en bestraffad-likelihood tillvägagångssätt, genom att bestraffa den totala variationen (TV) av logodds över grafen G. Detta är också känt som grafen smält lasso, eller noll-order graf trendfiltrering (Wang et al., 2014). Låt β (s) = log{w (s) /(1 − w (s) )} betecknar logodds som motsvarar vänster-split sannolikhet på plats s. Vi uppskattar {β (s) : s  till V } genom att lösa följande optimeringsproblem: där λ är en straffparameter vald adaptivt (se avsnitt 3.2). Detta är också känt som graf-fused lasso (GFL) under en binomial sannolikhet. För att upprepa: detta är problemet för en trädknuta. Vi löser ett separat sådant problem för varje icke-terminal nod, parallellt. Valet av 1 strafffunktion kan motiveras av både beräknings- och statistikskäl. Först leder det till ett problem med konvex optimering. De beräkningsmässiga fördelarna med konvexitet är en viktig faktor, eftersom problem (2) måste lösas hundratals eller tusentals gånger (en gång för varje trädsplit) på en stor rumslig lattice. För det andra, 1 straff leder till en icke-linjär smidigare, kan anpassa sig till olika grader av jämnhet i olika rumsliga regioner. Detta kontrasterar mot den linjära utjämning som ges av en 2 straff (eller motsvarande, en Gaussian CAR tidigare). Denna anpassningsförmåga har viktiga konsekvenser för uppskattningsnoggrannheten, både teoretiskt och praktiskt. Till exempel Sadhanala et al. (2016) har fastställt minimax-nivåer för att uppskatta funktioner på reguljära rutnät under squared-terror förlust, vilket utökar arbetet Donoho och Johnstone (1998) på 1D-funktioner. De visade att även om den totala variationen denoiserar minimax 9-hastigheten, misslyckas alla linjära jämnare (inklusive 2 - baserade jämnare) med att göra det. Dessa resultat gäller inte direkt vårt problem, som varken innebär ett reguljärt rutnät eller en squared-terror förlust. Men de är suggestiva att 1 straff kan föredras framför 2 straff för rumslig utjämning, om vi inte tror att det underliggande rumsliga fältet är homogent jämnt. (I vår applikation finns det starka skäl att tro motsatsen, på grund av line-of-sight ocklusions som kan sprida gammastrålar från fysiskt närliggande källor.) Figur 3 ger viss intuition om ett leksaksproblem. De fyra översta panelerna visar binomial graf-fusionerad lasso kontra det bakre medelvärdet från en Gaussian inneboende autoregressiv modell på ett problem där den sanna signalen är slät. De fyra nedre panelerna visar samma jämförelse för ett problem där den sanna signalen har en diskontinuitet. För det smidiga problemet, har 1 MAP uppskattningen bara något sämre rot medel-kvadrat fel än CAR-modellen, medan för problemet med diskontinuitet, har det en mycket bättre RMSE. Våra experiment senare i dokumentet vidareutvecklar denna punkt mer i detalj: förlusten av att använda en 1 påföljd på en jämn signal är liten, jämfört med vinsterna med att använda den när signalen har regioner med varierande jämnhet (inklusive diskontinuiteter). Vi inser att konvexa straff kommer med sina egna begränsningar, i att de har en tendens att överkrympa stora signaler (i detta fall stora hopp över kanterna i grafen) i närvaro av gleshet. Se Fan och Li (2001) och Carvalho et al. (2010) för de ofta förekommande respektive Bayesianska perspektiven på denna fråga, liksom det tidigare dokumentet av Geman och Reynolds (1992) för en diskussion om samma fråga i bildbehandling. Att införliva icke konvexa påföljder i rumslig utjämning är ett aktivt forskningsområde, och vi beskriver ett sådant tillvägagångssätt i avsnitt 3.3. Men detta kräver MCMC, som tyvärr inte skalas lika bra som vår metod. För att lösa (2) på skalan använder vi en variant av en ny metod för att lösa det graffunderade lassoproblemet baserat på trail decompositions (Tansey och Scott, 2015). Denna algoritm utnyttjar ett grundläggande resultat i grafteori som säger att alla anslutna grafer med 2k udda-graders hörn kan brytas ner i en uppsättning k icke-överlappande spår. Detta gör det möjligt för oss att skriva om straffterm i (2) som en summering över spår T: Ange den negativa logsannolikhet som (y, β), slack variabler z introduceras sedan för varje β (i) i straffet, vilket resulterar i följande motsvarande problem: 10 där begränsningarna innehar för alla par (r, s)  till τ, för alla τ  på T. Vi löser sedan (4) via alternativa riktningsmetod för multiplikatorer, eller ADMM (Boyd et al., 2011). Eftersom (y, β) inte är squared-terror förlust, en liten ändring av algoritmen i Tansey och Scott (2015) krävs. Som i Newtons metod, vid varje steg av algoritmen vi iterativt bildar en andra-order Taylor expansion till (y, β) till det aktuella värdet av β, och använda detta som en surrogat förlust funktion  på (., β). Detta ger följande uppdatering steg: där z (τ) är uppsättningen av slack variabler i τ th leden, u är den skalade ADMM dual variable, α är den skalar ADMM steg storlek parameter, J (s) är den uppsättning av slack-variable index som karta till β (s), och A är en gles binär matris som används för att koda lämpliga β (s) för varje z (j). Vi itererar uppdateringarna i (5)-(10) tills de dubbla och primära resterna har tillräckligt små normer. Observera att uppdateringen i (9) kräver att man löser ett vägt 1-dimensionellt smält lasso problem. Detta kan göras i linjär tid via en effektiv dynamisk programmeringsrutin (Johnson, 2013; Arnold et al., 2014). Vi väljer λ-hyperparametern separat för varje icke-terminal nod av dyadisk densitet nedbrytning, eftersom den nivå av utjämning som krävs på olika skalor kan skilja sig åt. Ett vanligt sätt att justera hyperparametrar är att använda korsvalidering. Det är dock inte klart i vårt fall hur vi ska dela upp våra data i utbildnings- och testdelningar så att korsvalidering skulle vara meningsfullt. Dessutom, eftersom vi använder en konvex optimering rutin, kan vi använda varma starter för att snabbt utforska lösningsvägen över ett rutnät av λ värden. Vi jämför sedan passformerna för olika λ's med hjälp av ett informationskriterium, att utnyttja resultaten på grad av frihet för 1 - baserad graf jämnare i Wang et al. (2014)................................................................ Vi använder BIC för vårt urvalskriterium, eftersom det är mer konservativt att hitta förändringspunkter än AIC eller AICc, även om resultaten är likartade i alla kriterier. I detta avsnitt beskriver vi en Bayesian variation av vår uppskattning. Detta tillvägagångssätt kan användas på problem av liten till måttlig storlek, och erbjuder tre fördelar jämfört med den klassiska 11 formulering: 1) Det möjliggör osäkerhet kvantifiering; 2) Det erbjuder en mer principfast metod för att välja hyperparametern λ; och 3) det erbjuder chansen att gå bortom 1 straff, via valet av tunga svans före som bevarar stora hopp i signalen över kanterna i grafen. Avyttringen är att det krävs MCMC, vilket leder till mycket längre beräkningstider. Vi observerar först att den graf-bearbetade lasso estimator som presenteras ovan är ett specialfall av en maximal a posteriori (MAP) uppskattning av β under den onormaliserade posterior där  och (K+1) är order-K graf trend filtrering straff matris (Wang et al., 2014). Konkret, när K = 0, matrisen  och (1) är den orienterade kantmatrisen för grafen G, som har en rad för varje kant e, och för vilken e,s = −1, och e,t = 0 för t /  till {s, r}. Denna matris extraherar den första ordningen skillnader i β över kanter, återhämta straff i (2). Högre ordning straff matriser (K > 1) definieras rekursivt, som i Wang et al. (2014)................................................................ I hela, vi använder d för att beteckna antalet rader i K+1. I (11), tidigare över β är en produkt av oberoende Laplace tidigare över (K + 1)-ordning skillnader i β över kanter. Denna tidigare kan representeras som en skala blandning av normals. Om vi hyr ut  på ett exponentiellt sätt, återvinns Laplace tidigare i (11) marginellt. Men om p( ska) är någon annan blandning distribution, vi får några andra global / lokal krympning innan på varje post i  och (K+1) β. Roualdes (2015) och Faulkner och Minin (2015) utforska detta ramverk för den endimensionella trend-filtrering fall, dvs.. där G är en kedja graf, förutsatt Gaussian fel. Vi utforskar samma klass av priors för binomial modeller på allmänna grafer, i samband med densitet utjämning problem. För fallet K = 0, som straffar första skillnader, det finns en nära koppling till Markov slumpmässiga fält, som den typ som används i Zhou och Hanson (2011). Eftersom varje diagonal post i  ska motsvara en kant (r, s) av grafen, kan vi skriva om (12) som, villkorat på λ och  ska, är ett Gaussian Markov slumpmässigt fält. Denna tidigare är felaktig på R p eftersom  och T är singular, men är korrekt på en subrymd definieras av en 12 lämplig begränsning. Specifikt, begränsa slutsatsen till subrymden underförstått av den centrerande begränsningen {β  på R p : 1 T β = 0}, där 1 är vektorn för alla 1's, säkerställer anständighet (c.f. Banerjee m.fl., 2004, Ch. Detta anslag är avsett att täcka utgifter för studier, möten med experter, möten med experter, möten med experter, möten med experter, möten med experter, möten med experter, möten med experter, möten med experter, möten med experter och experter samt andra administrativa utgifter. Detta följer av det faktum att A T är en viktad graf Laplacian matris, vars nullrymd spänns av 1 om alla diagonala poster i A är positiva. Analagiska begränsningar garanterar anständighet för högre ordning straff matriser, när K > 1. För mixningen distribution p(en) använder vi valet som leder till den generaliserade dubbel Pareto tidigare i Armagan et al. (2013), även kallad gamma-lasso tidigare beskriven i Taddy (2013). Detta val leder till tunga svansar och därmed undvika några av de overshrinkage problem av 1 /Laplace föregående. I litteraturen har man föreslagit många så tungstjärtade shrigage priors; se Polson och Scott (2011) för en genomgång. Vi använder gammalasso eftersom det också har en enkel tolkning som en blandning av Laplace tidigare, vilket gör kopplingen till grafen smält lasso explicit, och eftersom våra experiment med andra val ledde till i huvudsak likvärdiga prestanda. Vi representerar binomialsannolikheten med hjälp av Pólya-Gamma data förstärkningsschema som beskrivs i Polson et al. (2013)...................................................................... Detta ger villkorligt konjugera MCMC-uppdateringar för alla modellparametrar. För att upprätthålla anständighetsbegränsningen lägger vi i en global snap term till binomialprovtagningsmodellen, d.v.s. Avlyssna α tilldelas en vag tidigare, och vi centrerar β (s) är vid varje MCMC iteration. Den naturliga full Bayes-metoden för montering av denna modell är att inkludera den globala krympparametern λ 2 i provtagaren. Genom att göra det stötte vi på allvarliga problem, både med dålig blandning (en allmänt känd fråga med global / lokal krympning tidigare) och med numerisk stabilitet av matrisoperationer. Vi kunde inte få något av de MCMC system som beskrivs för 1D trend-filtrering modell av Roualdes (2015) och Faulkner och Minin (2015)-Gibbs provtagning och Hamiltonian Monte Carlo, respektive att ge stabil prestanda för allmänna grafer. Dessa frågor behandlades genom att montera separata modeller över ett diskret rutnät av fasta värden för λ 2, och välja λ 2 med hjälp av avvikelse informationskriteriet (Spiegelhalter et al., 2002). Vår strategi är alltså empiriska Bayes, snarare än helt Bayes. Fullständiga uppgifter finns i bilagan. Simuleringar Vi genomförde en rad experiment för att fastställa resultatet av vårt föreslagna tillvägagångssätt, inklusive scenarier där modellen är felspecificerad. I varje experiment använde vi ett 50 × 50 rutnät, där densiteten f (s) i varje rutnätscell var en Gaussisk distribution med 13 enheter varians och ett rumsligt varierande medelvärde μ (s) i intervallet [−3, 3]. Ett sådant problem är tillräckligt litet för att Bayesianska strategier ska kunna användas som riktmärken. Vi övervägde tre rumsliga beroendescenarier för medelvärdet baserat på bitvis polynom, för vilka vektorn (K+1) μ hade mestadels noll poster: bitvis konstant (K = 0), linjär (K = 1) och kvadratisk (K = 2). Dessa är mycket nästan splines av ordning 0 till 2, respektive; se Wang et al. (2014)................................................................ Vi använde också en fjärde, smidig signal med kontinuerliga derivat av alla beställningar, som en Gaussian CAR-modell bör passa. Figur 4 (vänster kolumn) visar de fyra olika grund sanningar för medelvärdet av varje Gaussian i våra experiment. Observationerna i varje rutnätscell diskretiserades till 2048 bins i intervallet [−6, 6]. Vi varierade antalet observationer per cell och mätte prestandan för varje metod, över alla rumsliga celler. För den empiriska Bayes-metoden använde vi en K = 2 trendfiltreringsmodell. Vi passar även Bayesian CAR modeller av MCMC, efter Zhou och Hanson (2011). För att mäta återuppbyggnadens noggrannhet beräknade vi det maximala CDF-felet mellan den verkliga fördelningen på plats s och den uppskattade fördelningen: där F och F är den sanna och uppskattade CDF, respektive. Skälet till denna felåtgärd är att den slutliga användningen för uppskattningen är i ett Kolmogorov-Smirnov test för anomalier, för vilka vi behöver CDF (se ansökan i avsnitt 5). Figur 4 (mellan- och högerkolonner) visar resultaten av dessa experiment. För bitvis konstant och bitvis linjära signaler, den binomial graf-smält lasso (GFL) leder till den bästa genomsnittliga densitet rekonstruktion, slå båda Bayesian modeller (andra kolumnen i figur). För de bitvis kvadratiska och släta signalerna har den binomiala GFL och Gaussian CAR-modellen jämförbar genomsnittlig prestanda, medan empirisk-Bayes trendfiltrering leder till bästa uppskattning. Den tredje kolumnen i figur 4 visar också att för alla tre bitvis polynomiska signaler har den GFL-baserade metoden det lägsta värsta felet (där "sämsta fallet" betyder en metods högsta CDF-fel i alla rumsliga celler). För den släta signalen, GFL-metoden har jämförbara värsta fall fel till CAR-modellen, och Bayesian graf trend-filtrering med hjälp av gamma-lasso tidigare har det bästa värsta-fall felet. Denna uppenbara nära-minimaxity av GFL-metoden är viktig för tillämpningen i avsnitt 5, där ett dåligt värsta fall fel kan leda till uppblåst falsk-larmfrekvens för radiologisk anomali upptäckt. Vi kommer också att se över denna punkt i avsnitt 4.2. Sammantaget verkar den empiriska Bayes-metoden vara starkast när den underliggande variationen är homogen, medan vårt föreslagna tillvägagångssätt baserat på den grafbearbetade lasso ger prestanda som är lika eller bättre för en CAR-modell, men med mycket bättre skalbarhet. Eftersom de empiriska Bayes- och CAR-modellerna kräver betydligt mer beräkningstid är det inte troligt att de är praktiska för skalning till stora datauppsättningar såsom radiologiska kartor i full stad. Till exempel, i dessa syntetiska experiment, uppskatta den rumsliga variationen för en enda split nod i binärt träd tar 23,4 sekunder för CAR-modellen i genomsnitt (för 10 000 14 MCMC prover), jämfört med 340 millisekunder för GFL algoritm-nära två storleksordningar snabbare. Dessutom måste varje algoritm köras hundratals gånger, en gång för varje delning nod i trädet. Följaktligen är den praktiska skillnaden mellan metoderna en av att få en densitetsberäkning för alla 2500 rumsliga celler i minuter, jämfört med timmar. Denna skillnad vidgas endast när grafens storlek ökar utöver 50 × 50 fall, och uppgiften att verifiera god MCMC blandning för CAR- och Bayesian GTF-modeller blir svårare. Vi genomförde ett andra simuleringsexperiment för att demonstrera skabbigheten hos vår metod, samt för att undersöka dess fördelar och nackdelar jämfört med andra skalbara metoder (särskilt Gaussian kärna utjämning). Som en mark-sanning scenario, skapade vi en 500 × 500 rutnät av densitet funktioner f (s), som är tillräckligt stor för att fånga finskalig struktur över centrum av en stad. Vi utformade detta scenario med två kriterier i åtanke: 1. Densiteterna f (s) bör vara nästan stationära (dvs. alla samma), så att kärnan utjämning bör leda till utmärkt prestanda för de flesta rumsliga celler. Detta utgör en utmaning för den metod för utjämning av densiteten som vi har föreslagit, och som inte är utformad för stationära scenarier. 2. Utgångspunkten är följande: Det bör finnas en rumsligt isolerad diskontinuitet, av det slag som kan orsakas av en line-of-sight ocklusion eller någon annan rumsligt lokaliserad radiologisk egenskap i den byggda miljön. Denna funktion ger en lokaliserad utmaning till Gaussian kärna utjämning, som antar smidig variation. För att uppfylla dessa kriterier konstruerade vi densitetsfunktionerna f (s) som en blandning av ett globalt bidrag f 0 plus platsspecifika avvikelser på grund av en isolerad radiologisk källa. Den globala bakgrundstätheten f0 valdes till det empiriskt observerade genomsnittliga bakgrundsspektrumet gamma-energi på alla platser på Pickle Research Campus, baserat på 18 timmars datainsamling. För att skapa små avvikelser från stationaritet injicerade vi sedan artificiellt två radioaktiva källor i rutnätet, i cellen (101,101): (1) en 100-milliCurie källa kobolt-60, och (2) en 100-milliCurie källa till cesium-137. De spektraldensiteter av den globala bakgrunden, kobolt-60, och cesium-137 visas i figur 5. Observera att dessa cesium- och koboltdensiteter (betecknade f Co och f Cs, respektive) inte är de kända gamma-ray emissionsspektra av dessa två isotoper. Snarare är de empiriska spektra uppskattas i en separat uppsättning experiment med samma detektor som används i våra bakgrundskollektiv experiment. Skillnaderna mellan det teoretiska spektrumet och det empiriska spektrumet av en radioaktiv isotop kan förklaras av många fysiska fenomen, såsom Comptonspridning. Eftersom vårt mål är att uppskatta vad en verklig data-insamling protokoll sannolikt skulle stöta på, använder vi den empiriska spektra snarare än den teoretiska spektra (se bilagan). Fotoner på varje plats antogs ha energi dras slumpmässigt från lokal densitet f (s), och att upptäckas med en genomsnittlig Poisson hastighet av λ (s) Cs. Varje term definieras enligt följande: • λ 0 är den globala bakgrundsfrekvensen, som antas vara 40 fotoner per sekund i alla celler, den genomsnittliga talhastighet som observerats under vårt fältexperiment. • λ Co är den genomsnittliga ankomsthastigheten för fotoner från koboltkällan på platsen s. • λ Cs är den genomsnittliga ankomsthastigheten för fotoner från cesiumkällan på platsen s. Co och λ Co på varje plats, använde vi följande ekvation: där mCi är storleken på källan i milliCuries, och d s är avståndet i meter från centrum av rutnätscell s till källan. Denna ekvation är härledd från välförståelig fysik, tillsammans med det fältexperiment som beskrivs i tillägget. Mark-truth lokala densiteter specificerades sedan via lämplig plats-specifik konvex kombination av bakgrund, kobolt och cesiumspektrum: där den gemensamma nämnaren Cs. Undantaget är att för alla platser nordväst om cellen (101,101) i rutnätet simulerade vi en line-of-sight ocklusion, vilket innebär att cesium och kobolt fotoner inte når dessa platser (dvs. är utspridda innan de kommer dit). Således f (s) (x) = f 0 (x) för dessa celler. Detta skapar en diskontinuitet nära cellen (101,101). Det är svårt att visualisera alla 250.000 densitetsfunktioner på en gång, men figur 6 visar dem för en delmängd av 64 celler: var tionde rad och kolumn i rutnätet från cellen (71,71) till cellen (141.141). Cellen (101,101) innehåller den isolerade strålningskällan. Det är lätt att se det rumsliga förfallet i bidraget från denna isolerade källa till spektra av närliggande celler, eftersom cesium har en stor topp i sitt spektrum runt energikanal 550. Det är också lätt att se effekten av ocklusionen, i de 9 mest nordvästra cellerna i figuren. Vi betonar att figur 6 zoomar in på den enda regionen i 500×500 rutnät där det finns någon rumslig variation. Effekten av den isolerade källan sönderfaller ganska snabbt när man rör sig bort från denna källa; vid cell (141.141) har cesiumtoppen nästan försvunnit. Överallt annars är densiteterna nästan stationära. Detta ger ett mycket gynnsamt scenario för Gaussian kärna utjämning (eftersom linjär krympning kommer att vara optimal över nästan hela grafen), och en utmaning för vår föreslagna metod. Vi simulerade 4 olika datamängder som spänner över en rad uppehållstider från 10 sekunder till 5 minuter. Med uppehållstid menar vi den teoretiska tid som en observatör har tillbringat med att samla in data i varje cell, med hjälp av ett instrument som det som används för att samla in data på University of Texas campus. Som med de verkliga data, vi bined alla observationer i 2048 jämnt fördelade energikanaler. Alla metoder fungerar på räkningarna i behållarna, snarare än med rå simulerade data. Med hjälp av dessa simulerade data jämförde vi den rumsliga densiteten i flera skalor mot två alternativ. Haar-Fisz-utjämning. Vi tillämpade den våget-baserade metoden av Fryzlewicz och Nason (2004) för att uppskatta densiteten oberoende i varje rumslig cell. Detta innebär ingen rumslig utjämning, men denoise observationer inom varje cell på ett sätt som försöker bevara både spisky och släta funktioner. Tröskelparametern för Haar-Fisz-metoden valdes av den giriga trädbaserade algoritmen som implementerades i haarfisz R-paketet (Fryzlewicz, 2009). Se Jansen (2006) för vidare diskussion och detaljer om metoden. Multiskale Gaussian kärna utjämning. Vi tillämpade också standardtekniken av Gaussian kärna utjämning för att ge en rumsligt jämn uppsättning uppskattningar för uppdelning sannolikheter (ekvation 1) i varje nod i multiscale träd nedbrytning. Gaussian kärna utjämning är i utbredd användning för skalbar rumslig utjämning i en mängd olika fält, från fMRI till klimatmodellering. I vårt fall är uppskattningen vid varje nod ersatt med ett lokalt vägt genomsnitt av sitt område, med varje grannes vikt proportionell mot exp(− 1 2 δ(s, t)/c) · nt, där s och t är nod platser som motsvarar målet och grannnoder, respektive δ är Euclidean avstånd funktion, n t är antalet försök som observerats vid nod t, och c är en bandbredd. Vi utforskade bandbreddsinställningar på 1, 1.5, 3 och 5 celler, vilket trunkerade kärnvikterna till noll efter 3c celler bort från målnoden för att ge en gles utjämningsmatris. Tyvärr, på ett problem detta stora (250 000 rumsliga platser), Bayesian metoder som tidigare beskrivs kan inte möjligt att tillämpas som riktmärken: de tar timmar att köra på ett problem med endast 2.500 noder, och deras beräkningskrav skala mycket snabbare än linjärt när antalet noder ökar. Figur 7 visar det genomsnittliga och värsta felet på alla platser för alla metoder, mätt i det tidigare försöket. Som väntat, eftersom sanningen är nästan stationär, Gaussian kärnan-smoothing uppskattningar har den lägsta genomsnittliga fel; högre bandbredder leder till mer utjämning, och presterar bättre i genomsnitt (vänster panel). Den högra panelen i figur 7 visar dock att det finns ett högt pris att betala för detta överlägsna genomsnitt 17 prestanda. Det värsta felet (dvs. Den värsta uppskattningen av alla rumsliga platser) i Gaussian-kärnan förblir nästan konstant även när data ackumuleras, oavsett kärnbandbredd. Dessutom, för Gaussiska metoder, det finns ett omvänt samband mellan genomsnittlig prestanda och värsta fall prestanda. Till exempel leder en bandbredd på 5 till den bästa genomsnittliga prestandan för alla celler. Men även med 5 minuters data, dess värsta rekonstruktion har ett fel på nästan 0,25, vilket är mycket dåligt med tanke på att den största möjliga CDF fel är 1. Det värsta felet i den GFL-baserade metoden, å andra sidan, blir mindre med mer data. Återvänder till Figur 6, Vi noterar att denna dåliga värsta-fall prestanda av Gaussian metoden inte är isolerad till en eller två celler. Panelerna B-D visar återuppbyggnadsfelet för varje metod för alla celler från (71,71) till (141, 141), i den region där det finns diskontinuitet. Både Haar-Fisz utjämning och grafen smält lasso metod har fel som är nästan rumsligt invariant i denna region. Den Gaussiska kärnan smidigare, å andra sidan, har mycket dålig prestanda längs hela kammen av diskontinuitet. Denna sämsta-fall prestanda blir inte bättre eftersom data ackumuleras. Figur 8 belyser ytterligare effekterna av avbrottet genom att zooma in på rekonstruktionerna vid två närliggande celler på motsatta sidor av ocklusionen, under en uppehållstid på 10s. Haar-Fisz modellen gör ingen utjämning alls här, eftersom dess adaptivt valda tröskelparameter behöll alla wavelet koefficienter i Haar-Fisz expansion av den underliggande densiteten. Den Gaussiska kärnan jämnare uppskattar en falsk spik i mitten av densiteten på (100,100), eftersom den inte kan anpassa sig till ocklusionen. GFL-metoden, å andra sidan, rekonstruerar båda densiteterna med hög trohet. På grundval av dessa experiment drar vi slutsatsen att GFL-metoden helt klart är det bästa valet för rumslig utjämning i storskalig radiologisk undersökning. Även i en (mycket ogynnsam) nära stationär miljö, är dess genomsnittliga rekonstruktion jämförbar med Gaussian kärna utjämning, om än något värre. Men dess värsta återuppbyggnadsfel är mycket bättre, på grund av oklarheten. Detta gör GFL-metoden mycket mer lämplig för användning i vår applikation för rumsligt medveten radiologisk anomalidetektering som beskrivs i nästa avsnitt. Global linjär krympning är helt klart olämpligt för detta problem: dess värsta fel kan inte åtgärdas med mer data (utom genom att inte utjämna alls), och det skapar därför en situation där falska alarm praktiskt taget garanteras i någon rumslig region med en kraftig förändring i bakgrunden. Radiologisk undersökning och anomalidetektering vid UT-Austin Som vi diskuterade i inledningen är huvudmotiveringen bakom vår flerskaliga rumsliga densitet smidigare att möjliggöra förbättrad upptäckt av radiologiska anomalier. Vi visar nu den föreslagna metoden uppnår detta mål med hjälp av de 18 timmars data från Pickle Research Campus som beskrivs i inledningen. Verkliga radiologiska avvikelser är sällsynta. Detta är bra ur folkhälso- och säkerhetssynpunkt. Men det utgör en allvarlig utmaning när man utformar en studie för att bedöma effektiviteten hos en metod för upptäckt av radiologiska anomalier. Utan budget eller personal för en stor fältstudie med faktiska radioaktiva material och komplicerade observationsprotokoll, är vi tvungna att förlita oss på simulerade anomalier. Vårt fokus ligger därför på att göra dem så realistiska som möjligt, genom att utnyttja två viktiga datakällor: 1) de 18 timmarna av normala bakgrundsdata för UT Pickle Research Campus, som beskrivs i inledningen, och 2) de experimentella data som samlas in på empiriska spektra av cesium och kobolt källor, som beskrivs i bilagan. Träning. Vi delade först upp campus bakgrundsdata i utbildning (-80 %) och tester (-20 %) uppsättningar. I synnerhet hade träningsdata 2 015 515 totala fotonantal, och testdata 546 064 totala fotonantal, över alla energikanaler. Vi delade upp campus i ett regelbundet rutnät på 50m × 50m rumsliga platser. Med hjälp av endast träningsdata konstruerade vi platsspecifika uppskattningar av gamma-ray-spektrumet via flerskalig täthetsutjämning, som beskrivs i avsnitten 2 och 3. Vi använde också träningsdata för att uppskatta en global (spatiellt invariant) bakgrundstäthet, helt enkelt genom att beräkna antalet i varje bin i alla träningsobservationer. Vi hänvisar till dessa som lokala respektive globala uppskattningar. Uppgifterna var långt ifrån jämnt fördelade över de rumsliga platserna: i utbildningsuppgifterna hade den genomsnittliga platsen 12,140 fotonräkningar, den uppgiftsrikaste platsen hade 249000 räkningar, medan den uppgiftsfattigaste platsen endast hade 64 räkningar. Figur 1 visar två exempel på uppskattningar för platser som är åtskilda med flera hundra meter, tillsammans med den globala uppskattningen och rådata på dessa platser. Utbildningsdata för dessa två platser uppvisar tydligt olika bakgrunder nära den breda toppen i spektrumet mellan kanalerna 50 och 200; jämför varje histogram med de vänstra två panelerna med hjälp av den streckade svarta linjen (den globala genomsnittliga bakgrunden) som en gemensam referenspunkt. Vårt förfarande rekonstruerar troget dessa skillnader utan den alltför stora variationen i det råa empiriska spektra. Jag testar. Antag att vi försöker simulera T sekunder av observation från plats s. Vi antar att, förutom den normala bakgrundsstrålning, gamma-ray spektrum också innehåller utsläppen från en onormal källa med gamma-ray spektrum f a, av storlek m milliCuries, och på avstånd d meter från observatören. För att generera en sådan observation använde vi testdata enligt följande. Ett viktigt inslag i denna experimentella design är att vi aldrig antar en känd bakgrund. (1) För att simulera bakgrundsbidraget, prov T en sekund bakgrundsobservationer genom att bootstrappa delmängd av testdata som kom från platsen s (dvs. provtagning av raderna med ersättning). Aggregera dessa observationer för att bilda en 19 enda simulerad T-sekund bidrag b (s) från bakgrunden:. Summan är komponentvis: både de bootstrappade observationerna x (s) [t] och summan är vektorer av längd 2048. (2) För att simulera anomaliens bidrag, använd ekvation (14) för att omvandla källstorleken s och avståndet d till en förväntad per sekund talhastighet λ. Ta en Poisson slumpmässig variabel n  på Poisson(Tλ) för T sekunder av observation. Prov n fotoner från den kända källtätheten f a (t.ex. de högra två panelerna i figur 5 ). Aggregera dessa prover i 2048 energikanaler för att bilda bidraget från den avvikande källan, a (s). (3) Kombinera bidragen från bakgrund och anomali, att bilda en enda "bootstrap + anomali" räkna vektor som simulerar T sekunder av observation på plats s. Målet med vår studie är att se hur enkelt vår metod kan skilja sådana anx (s) från normala (bakgrund-bara) observationer genom att jämföra det med bakgrundsspektrum som uppskattas från träningsdata. Riktmärken och simuleringsinställningar. Målet med detta dokument är inte att uppfinna en ny anomali-detektion algoritm, utan snarare att visa att rumslig utjämning kan förbättra befintliga algoritmer. Därför, för att avgöra om en simulerad (s) är en anomali, använder vi en ett-prov Kolmogorov-Smirnov test som jämför empirisk energi CDF avx (s) med CDF av den uppskattade bakgrundsspektrum på plats s. Detta är identiskt med metoden för Chan et al. (2014) med överspridningsparametern (= 1); generaliseringen till andra överspridningsparametrar skulle vara enkel i princip, även om vi inte utforskar detta. 2 För att fastställa att rumslig uppdelning och utjämning kan förbättra saker och ting, jämförde vi vår metod mot två andra: (1) Ett-prov KS/Chan et al. (2014) Test med hjälp av den övergripande uppskattningen av bakgrunden som referensfördelning. (2) KS-testet med två prover av Reinhart et al. (2015), där "prov 1" omfattar utbildningsobservationer från plats s, och "prov 2" är simulerad x (s). 2 Observera att författarna till Chan et al. (2014) förutsätter helt enkelt att bakgrunden är känd. De simuleringsparametrar som ska varieras är källstorleken m, källavståndet d, källspektrumet f a, och observationstiden T. Kartläggningen från källstorlek och avstånd till förväntad räkningshastighet (ekvation 14) är dock många till en. (Den översta panelen i figur 9 visar detta förhållande.) Vi parametriserar därför våra resultat i termer av den förväntade anomalien talfrekvensen snarare än storlek och avstånd individuellt. Vi simulerade data över alla kombinationer av följande inställningar: • T-tid {2, 4, 6, 10, 14, 20, 30, 40, 60, 90, 120, 180} sekunder. • anomaliräkningsfrekvens r 0, 1, 0, 25, 0, 5, 1, 2, 5, 10, 15, 20, 25, 50} fotoner per sekund. ............................................................................................................................... Se figur 5. För varje kombination av inställningar, använde vi bootstrap-baserade procedur bara beskrivs för att simulera 1000 anomalous observationerx (s) på varje plats. För ett intervall av detektionströsklar t beräknade vi den verkliga positiva hastigheten för varje förfarande, dvs. fraktionen av simulerade x (s) vars KS-statistik överskred t. Vi beräknade också den falska positiva hastigheten, dvs. Sannolikheten för att KS-statistiken för en normal observation kommer att överskrida t. För varje metod aggregerade vi de falska positiva och sanna positiva värdena på alla rumsliga platser för att beräkna en metodspecifik ROC-kurva för varje simuleringsinställning. Val av nollhypotes. För att vara försiktiga beräknade vi inte den felaktiga positiva frekvensen för varje metod med hjälp av teoretiska referensfördelningar för KS-statistiken för ett eller två stickprov. Snarare beräknade vi en "empirisk null" för varje metod genom att tillämpa varje KS-test på testdata utan artificiella anomalier injiceras. Dessa empiriska nullfördelningar, tillsammans med standard KS referensfördelning, visas i de fyra nedre panelerna i figur 9. De empiriska nullorna ligger nära den teoretiska nullen men har något tjockare högerstjärtar, eftersom KS-testet med ett enda prov inte ger upphov till osäkerhet i bakgrundsuppskattningarna. Genom alla simuleringsinställningar gav KS-testet baserat på rumsligt utjämnade densitetsberäkningar jämnt bättre ROC-kurvor än de andra metoderna. KS-testet med ett enda prov baserat på en uppskattning av den globala densiteten utför det värsta, medan KS-testet med två prov är mellanliggande mellan de två. Detta visar fint att det finns två tydliga källor till fördelar här: en prestanda bula på grund av rumslig upplösning, och en annan tydlig bula på grund av rumslig utjämning. Figur 10 ger ett litet urval av våra resultat. I dessa tre paneler är den förväntade räkningshastigheten från den anomaliska källan 5 gammastrålar per sekund, jämfört med en bakgrundshastighet på 39 per sekund, och observationstiden ökar från 20 till 40 till 60 sekunder avläsning vänster till höger. En noggrann jämförelse av panelerna visar att KS-testet med lokala 21 (spatiellt utjämnade) densitetsberäkningar i huvudsak fungerar lika bra med 40 sekunders data som metoden med två prover gör med 60 sekunders data, vilket tyder på en ungefär 30-procentig förbättring av tid till upptäckt. Liknande förbättringar sågs i alla simuleringsinställningar. Vi har presenterat multiscale rumslig densitet utjämning som en ny teknik för att uppskatta en rumsligt varierande densitet funktion. Våra bevis tyder på att metoden erbjuder toppmoderna prestanda när den införlivas i anomali-detektion protokoll för radiologiska undersökningsdata. De kraftförbättringar vi fann i förhållande till befintliga metoder är inga små praktiska frågor, särskilt i ett brottsbekämpande sammanhang (där upptäckt av en källa med högt självförtroende 30 % snabbare kan göra stor skillnad). Dessutom visar simuleringsbeviset i avsnitt 4 tydligt att, åtminstone för radiologiska undersökningar, rumslig utjämning med vår metod ger en gynnsam blandning av bra genomsnittliga fel med bra värsta vård fel. Det är också tydligt att global linjär krympning (antingen i form av en CAR-modell eller Gaussisk kärnutjämning) inte kan hantera den typ av skarpa förändringar i bakgrunden som förekommer i verkliga radiologiska data. En potentiell brist är att även om vår metod jämnar ut densiteten i rymden mycket effektivt, så jämnar den inte ut densiteten över dess underliggande stöd B på ett särskilt sofistikerat sätt. (Som det nu är, det enda sättet att införliva denna form av utjämning är att använda en grundare trädsammansättning, som har en effekt liknande den att använda en bredare bin bredd i ett histogram.) Ett självklart sätt att generalisera metoden skulle alltså vara att slå samman vårt arbete med några av de idéer som diskuteras av exempelvis Jansen (2006) och Willett och Nowak (2007), för att jämna ut energispektrumet. (Ingen av dessa tekniker inkluderar rumslig utjämning.) Det är inte omedelbart uppenbart för oss hur vi ska åstadkomma denna syntes, men detta är ett aktivt forskningsområde. Tackar. Författarna tackar Patrick Vetter från UT Applied Research Laboratories för hans hjälp med de pilotstudier som beskrivs här; University of Texas Police Department för deras pågående samarbete med datainsamling; och Ryan Tibshirani från CMU för att dela med sig av sin expertis om algoritmer för graf smält lasso.
För det första, för bakgrund om rumslig densitet utjämning över grafer, hänvisar vi läsaren till papperet genom REF, samt referenser däri.
Multiscale spatial density smoothing: an application to large-scale radiological survey and anomaly detection
{'venue': 'Journal of the American Statistical Association, vol. 112 no. 519 (2017), pp. 1047-1063', 'journal': None, 'mag_field_of_study': ['Mathematics', 'Physics']}