Eötvös Loránd Tudományegyetem Matematika Intézet
Daróczy-Kiss Endre
Galaxis-eloszlások vizsgálata BSc szakdolgozat Témavezető: Pröhle Tamás
ELTE Valószínűségelméleti és Statisztika Tanszék 2015. Budapest
Köszönetnyilvánítás Ezúton szeretnék köszönetet mondani Pröhle Tamás Professzor Úrnak a témavezetéséért, a segítőkészségéért, a türelméért, és hogy bármikor bekopogtam hozzá, mindig tudott időt szakítani a kérdéseim megválaszolására. Szeretném megköszönni Balázs Lajos Professzor Úrnak, hogy segített megtalálni ezt az izgalmas témát, és a sok segítséget a csillagászati fogalmak megértésében és ezen fejezetek leellenőrzésében. Köszönettel tartozom Rácz Istvánnak, aki a Millennium szimuláció működését mutatta meg. A legtöbb köszönet azonban az egész családomnak de különösképpen a szüleimnek és a nagyszüleimnek jár, a sok támogatásért, bíztatásért és szeretetért, ami az egyetemi éveimet vidámabbá és gondtalanabbá tette. Szeretném megköszönni Nórinak is a sok támogatást, bíztatást, és hogy lelket öntött belém a reménytelen pillanatokban is.
iii
Tartalomjegyzék 1. Kozmológiai háttér
2
1.1. Univerzumunk dióhéjban . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2. GRB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.3. A Millennium szimuláció . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.4. Az elvégzendő vizsgálatok . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2. Matematikai módszerek
7
2.1. Pontfolyamatok: Poisson-folyamat . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2. Térbeli véletlenszerűség . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.3. CSR tesztek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.4. Megjegyzések a fejezethez . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
3. Elemzések
24
3.1. A vizsgálni kívánt adatok leírása . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.2. Távolságon alapuló tesztek . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.3. A
χ2 -próba
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
4. Összefoglalás és további kérdések
28
Irodalomjegyzék
29
iv
Ábrák jegyzéke 1.1. Ősrobbanás [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2. GRB eloszlás [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.3. Sötét anyag eloszlás, nagy és kisebb léptékben [3], [4].
. . . . . . . . . . . . .
5
2.1. Beérkezési idők . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.2. Si val. változók . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.3. Az (a, b] intervallumba eső események száma . . . . . . . . . . . . . . . . . . .
10
2.4. Néhány térbeli eloszlás [10] . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.5. Határ probléma [21] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
2.6. Voronoi felbontás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.1. A G függvény, szimultán MC-vel . . . . . . . . . . . . . . . . . . . . . . . . .
25
3.2. Az F függvény, szimultán MC-vel . . . . . . . . . . . . . . . . . . . . . . . . .
26
3.3. A K függvény, szimultán MC-vel . . . . . . . . . . . . . . . . . . . . . . . . .
27
v
Bevezetés A modern fizikai kozmológiában, a kozmológiai elv szerint, kellően nagy léptékben vizsgálva a Világegyetemet azt tapasztaljuk, hogy az anyag eloszlása homogén és izotrop [16]. Ez felületesen megfogalmazva a következőt jelenti: az Univerzumban jelenlévő anyag átlagos sűrűsége mindenhol ugyanannyi, és bármilyen irányban vizsgáljuk a Világegyetemet, az ugyanolyannak mutatkozik. Ezt az axiómát látszik cáfolni, I. Horváth, J. Hakkila és Zs. Bagoly által nemrég publikált cikk [12], miszerint a Földtől több milliárd fényévre lévő objektumok csomósodást mutatnak. Ilyen nagy távolságban, speciális teleszkópok segítségével, csak olyan objektumokat tudunk megfigyelni, amelyek nagy energiát bocsájtanak ki. Ezek közül a rövid ideg tartó nagy energia kibocsátással járó jelenséget gamma kitörésnek, vagy GRB-nek (gamma-ray burst) nevezzük [18]. A gamma kitörések, túlnyomó részt nagytömegű csillagok halálakor jönnek létre. Nagy tömegű csillagok viszont általában olyan galaxisokban keletkeznek, amelyek csillagkeletkezésben aktívak. A dolgozatban a csillagkeletkezésben aktív galaxisok térbeli eloszlását vizsgálom, statisztikai módszerekkel, amihez az adatokat a Millennium szimuláció szolgáltatja. A dolgozat felépítése a következő: az 1. fejezetben ismeretterjesztői szinten leírom az Univerzum- és a GRB-nek a feladat szempontjából fontos tulajdonságait, néhány szót ejtek a Millennium szimulációról, és pontosítom az elvégzendő feladatot. A 2. fejezetben részletezem az alkalmazott matematikai módszereket, a 3. fejezetben pedig az elvégzem az elemzéseket az R programcsomag segítségével.
1
1. fejezet
Kozmológiai háttér 1.1. Univerzumunk dióhéjban Jelenleg a legszélesebb körben elfogadott modell, ami az Univerzumunk kezdetét és fejlődését leírja, az ősrobbanás elmélete [20], [15]. Eszerint a világegyetem egy rendkívül forró és sűrű állapotból fejlődött ki. Az ősrobbanás idővonalát szemlélve az Univerzum keletkezése egy végtelen hőmérsékletű és sűrűségű, nulla sugarú gömbbel, a szingularitással kezdődött. Ezt követte 0-tól 10−43 másodpercig a Planck-korszak, ahol az univerzumot hatalmas energiasűrűség és az ezzel együtt járó óriási nyomás és hőmérséklet töltötte ki egyenletesen és izotrop módon. A következő szakasz a kozmikus infláció, amikor a világegyetem exponenciális növekedésbe kezd. Ahogy az univerzum tágul, a hőmérséklete csökken. Nagyjából a 380 000. évben létrejönnek az első atommagok melynek következtében a Világegyetem átlátszóvá válik, megjelennek a szabad fotonok. Ezt a fénylést a mai napig meg lehet figyelni a kozmikus háttérsugárzás képében. Ennek a sugárzásnak a hőmérséklete az égbolton megközelítőleg egyenletes. Az idő múlásával a nagyjából egyenletes eloszlású anyag kissé sűrűbb részeinek gravitációs vonzása magukhoz húzta a környező anyagot, amelyek egyre sűrűbbé váltak, létrehozva ködöket, csillagokat, galaxisokat és egyéb csillagászati struktúrákat. Ennek a folyamatnak a részleteit az Univerzumban jelenlévő háromféle anyag, a hideg sötét anyag, a forró sötét anyag és a barionos anyag mennyisége határozza meg. A sötét anyag nem megfigyelhető, semmilyen elektromágneses sugárzást nem bocsájt ki és nem nyel el, jelenlétére csak az anyagra kifejtett gravitációs hatásból következtetünk. Legjobb tudásunk szerint az Univerzumban található anyag részben a sötét anyag gravitációs hatására sűrűsödött galaxisokká. Legpontosabb méréseink azt mutatják, hogy az univerzum tömegében a sötét anyag aránya 22%, a megfigyelhető anyag csupán 4% , a maradék 74% pedig sötét energia. A sötét energia az a feltételezett energia forma, amely az egész Univerzumban jelen van, és erős antigravitációs hatást fejt ki. Ez a legelfogadottabb elmélet annak magyarázatára, hogy a világegyetem 2
gyorsulva tágul. Mindezen elméleteknek az egyik matematikai modellje a ΛCDM (Lambda-Cold Dark Matter) modell.
1.1. ábra. Ősrobbanás [1]
1.2. GRB A gamma-kitörések felfedezése A gamma-sugár nagy frekvenciájú, rövid hullámhosszú, nagy energiájú elektromágneses sugárzás, ami szabad szemmel nem látható. A természetben is keletkezik, és mesterségesen is előállítható [17]. A gamma-kitörések (rövidítve : GRB, az angol gamma-ray burst-ből) felfedezése [7] a véletlen művének tulajdonítható. A hidegháború idején egy amerikai kutatócsoport feladata az atombomba-robbanások detektálása volt. Ezt, az ilyen robbanásokat követő gamma-sugarak világűrből való észlelésével próbálták figyelni. Sok gamma-felvillanást érzékeltek, de nem a Föld irányából, hanem a Világűrből.
3
A GRB-k helyeinek meghatározása sokáig problémát jelentett. A hetvenes évek végén, nyolcvanas évek elején felbocsátott szovjet Venera 11, 12, 13 és 14 holdakon elhelyezett Konus műszerek után az 1991-ben üzembe állított CGRO (Compton Gamma Ray Observatory) olyan műhold volt, amelyet kifejezetten a gammakitörések vizsgálatára hoztak létre. A több nagy látószögű detektorral felszerelt műhold már csaknem háromezer GRB helyét tudta meghatározni.
A gamma-kitörések keletkezése A GRB-k keletkezését egyelőre homály fedi, azonban komoly gyanúsítottak a legalább 40 naptömegű, saját tengelyük körül gyorsan forgó csillagok, melyek egy szupernova robbanásban fejezik be életüket [7], [11]. Más elméletek szerint kettős neutroncsillagok (nagy tömegű csillagok végállapotai, melyek nagy mennyiségű neutront tartalmaznak) összeütközése is lehet a gamma-kitörés kiváltó oka. Fontos megjegyezni, hogy a fent említett keletkezési módok túlnyomó többsége a nagytömegű csillagokból indul ki. Nagytömegű csillagok viszont egy galaxis csillagkeletkezési régiójában jönnek létre, tehát ahol ilyen gamma kitörést látunk, feltételezhetően ott egy galaxis csillagkeletkezésben aktív régiója is megtalálható.
1.2. ábra. GRB eloszlás [2]
A gamma-kitörések jelentősége Ahogy azt az előző alfejezetben olvashattuk, a GRB-k nagy tömegű csillagok halálakor jönnek létre, amiből következik, hogy a megfigyelt GRB egy galaxis közelében található. Így a GRB-k térbeli eloszlása segíthet az Univerzum nagyléptékű szerkezetének a megismerésében. A fent említett cikkből kiderül, hogy 283 GRB vöröseltolódását sikerült megmérni (a vörös eltolódás és a távolság között korreláció van). A megközelítőleg 300 megfigyelésből, a körülbelül 9 és 11 milliárd fényévre található 31 GRB eloszlása nem egyenletes. Közel a 4
felük 14-15 kitörés az égbolt egy nyolcadán található. A cikk azt is megmutatja, hogy egy ilyen csomósodás véletlen kialakulása nagyon valószínűtlen. Az 1.2 ábrán a kék pontok a 283 GRB-t, a piros pontok a nagy távolságban lévő 31 gamma-kitörést jelölik.
1.3. A Millennium szimuláció A Millennium szimuláció egy, a Max Planck Institute For Astrophysics által futatott kozmológiai N-test szimuláció (fizikai erők hatása alatt lévő részecskék dinamikai rendszerének szimulációja), amely 1010 nagyságrendű elemet tartalmaz, egy 2 milliárd fényév oldalú kockában [13], [19]. Célja a sötét anyag eloszlásának szimulálása, de nyomon követhetjük a körülbelül ∼ 107 nagyságrendű, a Kis Magellán Felhőnél fényesebb galaxisok fejlődését is. A keletkezett adatokat egy mindenki számára hozzáférhető adatbázisban tárolják. A szimuláció a ΛCDM matematikai modell szerint működik.
1.3. ábra. Sötét anyag eloszlás, nagy és kisebb léptékben [3], [4].
1.4. Az elvégzendő vizsgálatok A Millennium szimuláció adatbázisából kiválasztunk véletlenszerűen 104 nagyságrendű, csillagkeletkezésében aktív galaxist. Minden galaxishoz a következő paraméterek tartoznak: – galaxyID: ezzel a számmal azonosítjuk a különböző galaxisokat, – (x,y,z): a galaxis koordinátái, A galaxisok térbeli elhelyezkedését pontfolyamatnak tekintjük. Az előzőekben megfogalmazottak és a fent említett [12] cikk szerint a GRB-k eloszlása nem véletlenszerű. Ezért azt tűztük ki célunknak, hogy a GRB-kel feltételezhetően kapcsolatban álló, csillagkeletkezésben aktív galaxisok eloszlásáról egy állítást tudjunk kimondani.
5
2. fejezet
Matematikai módszerek 2.1. Pontfolyamatok : Poisson-folyamat A pontfolyamatok elmélete az alkalmazott statisztika egy olyan eszköze, amellyel térbeli adatokat lehet elemezni. Az asztrofizikában többek között ezt a jól kidolgozott módszert használják, például galaxis eloszlások vizsgálatára. Ebben az alfejezetben a legegyszerűbb pont folyamatot, a Poisson-folyamatot fogjuk áttekinteni.
2.1.1. Fontosabb eloszlások Mielőtt belekezdenénk a pontfolyamatok tárgyalásába, leírom a Poisson-folyamatban megjelenő fontosabb eloszlásokat. Binomiális eloszlás A binomiális eloszlású diszkrét valószínűségi változóval a visszatevéses mintavételt lehet modellezni, tehát, ha elvégzünk n-szer egy véletlen kísérletet, akkor megfigyeljük, hogy egy adott esemény hányszor következik be. Ennek az adott kísérletnek a valószínűsége: p, a bekövetkezett esemény gyakorisága : k(= 1,2, . . . n). 1. Definíció. Egy X valószínűségi változó n és p paraméterű binomiális eloszlást követ, akkor és csak akkor, ha eloszlás függvénye : n k P (X = k) = p (1 − p)n−k , k
k = 0, 1, 2, ..., n
– A binomiális eloszlás várható értéke : E(X) = np p – A binomiális eloszlás szórása: D(X) = np(1 − p) Poisson eloszlás A Poisson-eloszlás egy diszkrét valószínűségi eloszlás, ami adott idő alatt, az ismert valószínűséggel megtörténő események bekövetkezésének számát írja le. 7
2. Definíció. Egy X valószínűségi változó, β paraméterű Poisson-eloszlást követ, akkor és csak akkor, ha eloszlásfüggvénye : P (X = k) =
β k −β e , k!
k = 0, 1, 2, ...
alakú, ahol β > 0 konstans. – A Poisson-eloszlás várható értéke : E(X) = β √ – A Poisson-eloszlás szórása: D(X) = β Exponenciális eloszlás 3. Definíció. Az X valószínűségi változó β paraméterű exponenciális eloszlást követ, ha sűrűségfüggvénye: f (x) =
βe−βx
ha x ≥ 0,
0
ha x < 0
– Az exponenciális eloszlás várható értéke: E(X) = – Az exponenciális eloszlás szórása : D(X) =
1 β
1 β
– Az exponenciális eloszlás eloszlásfüggvénye: 1 − e−βx F (x) = 0
ha x ≥ 0, ha x < 0
Egyenletes eloszlás 4. Definíció. Az X valószínűségi változó egyenletes eloszlást követ az [a, b] intervallumon, ha sűrűségfüggvénye : f (x) =
1 b−a
ha a ≤ x ≤ b,
ha b < x < a
0
– Az egyenletes eloszlás várható értéke: E(X) = – Az egyenletes eloszlás szórása : D(X) =
a+b 2
b−a √ 12
– Az egyenletes eloszlás eloszlásfüggvénye: 0 F (x) = x−a b−a 1 8
ha x < a, ha a ≤ x ≤ b, ha x > b
Többdimenziós egyenletes eloszlás 5. Definíció. Az X valószínűségi változó a B ⊂ Rd (ahol d = 2 vagy 3) halmazon egyenletes eloszlást követ, ha sűrűségfüggvénye a következő : λd 1(G) ha x ∈ G f (x) = 0 ha x ∈ /G itt λd (G) a G halmaz területét, vagy térfogatát jelöli. Gamma eloszlás 6. Definíció. Az X valószínűségi változó α rendű, β paraméterű gamma-eloszlást követ, akkor és csak akkor, ha a sűrűségfüggvénye : α β xα−1 e−βx Γ(α) f (x) = 0 R ∞ α−1 −t e dt a gamma függvény. ahol Γ(α) = 0 t
ha x > 0, ha x ≤ 0
– A gamma-eloszlás várható értéke: E(X) = αβ , √
– A gamma-eloszlás szórása : D(X) =
α β ,
– Ha α = 1 akkor a gamma-eloszlás megegyezik a β paraméterű exponenciális eloszlással, – Ha α =
n 2
és β = 12 , akkor a gamma eloszlás megegyezik az n szabadságfokú χ2 elosz-
lással
2.1.2. Szemléletes bevezetés a pontfolyamatokba Az egy dimenziós pontfolyamat jól modellezi a véletlenszerű időpontokban bekövetkezett események egy sorozatát, például a véletlenszerű időpontokban beérkezett telefonhívásokat egy telefonközpontba. Ha egy térképen jelezzük azoknak az embereknek a helyzetét, akik a telefonközpontba telefonáltak, akkor egy olyan kétdimenziós pontmintához jutunk, ahol a telefonálók helyzete véletlenszerű. Az ilyen véletlen pontmintákat jól lehet modellezni ddimenziós pontfolyamatként, ahol d ≥ 2. T1
T2
T3
T4
2.1. ábra. Beérkezési idők Az egydimenziós pontfolyamatot matematikailag sokféleképpen lehet kezelni. Az egyik ilyen direkt módszer, ha a T1 < T2 < T3 . . . beérkezési időket vizsgáljuk, ahol Ti olyan valószínűségi változó, ami azt az időpontot jelöli, amikor az i. esemény (telefonhívás) bekövetkezik [8]. 9
Másik lehetőség ha a pontfolyamatot az Si val. változóval jellemezzük, ahol Si = Ti+1 −Ti . Az Si használata azért kényelmes ellentétben a Ti -vel mert bizonyos modellekben (pl.: Poisson folyamat), S1 , S2 , S3 . . . val. változók függetlenek, míg a Ti val. változók nem függetlenek a Ti < Ti+1 tulajdonság miatt. S1
S2
S3
S4
S5
2.2. ábra. Si val. változók A pontfolyamatot az Nt val. változó is jól leírja, ahol Nt = a t időpontig bekövetkezett események száma, azaz Nt =
∞ X
1{Ti ≤ t}
i=1
minden t ≥ 0, ahol 1{. . .} egyenlő 1-el ha a {. . .} zárójelek közötti állítás igaz, máskülönben 0-val egyenlő. Utoljára, nézzük az N (a, b] = Nb − Na val. változót. Itt N (a, b] (0 ≤ a ≤ b) az (a, b] intervallumba eső események számát adja meg. Poisson folyamatnál ez a módszer diszjunkt intervallumokra független.
a
b N (a, b] = 2
2.3. ábra. Az (a, b] intervallumba eső események száma Magasabb dimenziókban nincsen rendezés, ezért nincsen természetes analógiája az Si vagy az Nt val. változóknak. Viszont N (a, b]-nek van, ezért a legkézenfekvőbb ezt a fajta jellemzését a pontfolyamatoknak kiterjeszteni [8]. Legyen B ⊂ Rd korlátos, zárt halmaz, és x1 , x2 , . . . , xn ∈ Rd pedig a pontfolyamat pontjai. Ekkor N (B) =
n X
1{xi ∈ B}
i=1
azaz, N (B) a B halmazba eső pontok számát adja meg. Néha elegendő a V (B) indikátort használni a pontfolyamat jellemzésére, ahol V (B) = 1{N (B) = 0}
10
2.1.3. Egydimenziós Poisson-folyamat 7. Definíció. Egyenletes β > 0 intenzitású, egydimenziós R-beli Poisson-folyamatnak nevezzük azt a pontfolyamatot, ami az alábbi feltételeket kielégíti [8], [10]: (T1) minden korlátos, (a, b] intervallumra, az N (a, b] val. változó Poisson eloszlású β(b − a) paraméterrel, azaz P (N [a, b) = k) =
(β(b − a))k −β(b−a) e k!
k = 0,1, . . . ,
(T2) ha az (a1 , b1 ), . . . , (am , bm ) korlátos intervallumok diszjunktak, akkor az N (a1 , b1 ], . . . , N (am , bm ] val. változók függetlenek. Érdemes megjegyezni még néhány fontos tulajdonságát az R-beli Poisson-folyamatnak : 1. Az Si val. változók exponenciális eloszlásúak, β paraméterrel: P (Si ≤ s) = 1 − e−βs 2. Az Si val. változók függetlenek. 3. A Ti val. változók gamma eloszlást követnek α = i és β paraméterrel: f (t) =
β α α−1 −βt t e Γ(α)
minden t > 0 számra,és minden más esetben nullával egyenlő, 4. Nem történnek szimultán események, azaz P (N (a, a + h] ≥ 2) = o(h) ha h & 0, 5. Egy adott (a, b] intervallumba eső pontok számának a várható érteke arányos az intervallum hosszával: EN (a, b] = β(b − a), ahol β > 0 a folyamat intenzitása. 6. Bármely (a, b] korlátos intervallumban az események eloszlása egyenletes. Megjegyzés A Poisson-folyamatnak a fenti esetnél egy általánosabb formája az amikor a β intenzitás időben változik [8]. Ekkor: Z EN (a, b] =
b
β(t) dt a
ahol β(t) > 0 az intenzitás függvény. Ezt a folyamatot inhomogén Poisson folyamatnak nevezzük, és itt is igaz az a tulajdonság, hogy diszjunkt intervallumokra az N (a, b] függetlenek. 11
2.1.4. Térbeli Poisson-pontfolyamat 8. Definíció. Az Rd -beli (esetünkben d = 2 vagy 3) pontfolyamat, egyenletes, β > 0 intenzitású térbeli Poisson-pontfolyamat, ha a következő feltételeket teljesíti [8] : (TT1) minden korlátos, zárt B halmazra, az N (B) val. változónak Poisson-eloszlása van βλd (B) paraméterrel, (TT2) ha a B1 , B2 , . . . , Bm diszjunkt halmazok, akkor az N (B1 ), N (B2 ), . . . , N (Bm ) val. változók függetlenek. Itt a λd (B) a B halmaz területét illetve térfogatát jelöli. A fenti (TT1) és (TT2) feltétel egyértelműen meghatározza a térbeli Poisson folyamatot. Az egydimenziós esethez hasonlóan a többdimenziós Poisson folyamat néhány fontos tulajdonsága: 1. EN (B) = βλd (B) 2. P (N (B) > 1) = o(λd (B)), minden kicsi λd (B)-re 3. Legyen N (B) = n, ekkor az n darab esemény független véletlen minta a többdimenziós egyenletes eloszlásból B-n.
12
2.2. Térbeli véletlenszerűség Ebben az alfejezetben azt a kérdést próbáljuk megválaszolni, hogy hogyan néznek ki azok a pontminták, amelyek véletlenszerűen oszlanak el a térben.
2.2.1. Teljes térbeli véletlenszerűség Számos gyakorlati alkalmazásban találkozunk olyan feladattal, amelyeket térbeli ponthalmazként lehet modellezni, ahol a pontokat eseményekkel azonosítjuk. Ilyenek például: egy erdőben az egyes fák helyzetének a vizsgálata, galaxisok tömegeloszlásának modellezése, geológiai adatok elemzése. Sokszor az ilyen típusú feladatoknál fontos annak az eldöntése, hogy a kapott pontminta eloszlása véletlenszerű-e ? Vizsgáljuk meg az alábbi három ábrát, amelyeket közvetlenül el tudnunk érni az R programban. Az 1. ábrán, japán feketefenyők eloszlását láthatjuk, egy 5.7 méter oldalú négyzetben, a 2. ábrán vörösfenyők eloszlását láthatjuk egy 23 méter oldalú négyzetben, a 3. ábrán pedig sejtek eloszlását az egységnégyzeten. Intuíciónk azt súgallja, hogy az 1. ponteloszlás véletlenszerű. J. Diggle nyomán, a 2. ponteloszlást klaszterezett, a 3.-at pedig pedig szabályos pontmintának nevezzük [10]. A kérdés az, hogy a véletlenszerű ponteloszlást hogyan lehetne formalizálni ? Ezt vizsgáljuk a következőkben.
2.4. ábra. Néhány térbeli eloszlás [10] Ha egy pontminta véletlenszerűen oszlik el a térben akkor a következő tulajdonságok teljesülnek rá [10] [8]: 1. Az események száma egy B ⊂ Rd (d = 2 vagy 3) régióban ahol λd (B) a d dimenziós térfogatot jelöli, Poisson eloszlást követ βλd (B) paraméterrel, ahol a β konstans az intenzitást jelöli. 2. Az n darab esemény a B régión független minta a d-dimenziós egyenletes eloszlásból. A fenti megállapítások miatt a teljes térbeli véletlenszerűség (complete spatial randomness, rövidítve: CSR, és innentől kezdve CSR-ként hivatkozunk rá) matematikai modellje a térbeli homogén Poisson folyamat. Ahhoz, hogy a CSR hipotézisünket vizsgálhassuk, teszteket kell elvégeznünk. A következőkben ezeket a teszteket tárgyaljuk. 13
2.3. CSR tesztek Ebben a fejezetben a teljes térbeli véletlenszerűségre, azaz CSR-re nézünk át néhány olyan tesztet, amelyek a pontmintában lévő pontok távolságain alapulnak. Megnézünk ezen kívül egy olyan tesztet is, ami a távolságmódszerek egy alternatívája. Mindezek előtt azonban átnézzünk három Monte Carlo tesztet, és egy becslést a homogén Poisson-folyamat β intenzitására.
2.3.1. A β paraméter becslése A CSR tesztek menete általában a következőképpen folyik: Megfogalmazzuk a nullhipotézisünket miszerint a pontminták CSR, feltesszük, hogy H0 igaz, veszünk egy V teszt statisztikát, meghatározzuk a Vˆ tapasztalati eloszlásfüggvényt, és összehasonlítjuk az elméleti eloszlásfüggvénnyel. Ha az eltérés a két függvény között "nagy", akkor H0 -t elvetjük. Ezt a módszert a következő alfejezetekben részletesebben tárgyaljuk. Az adott elméleti eloszlásfüggvényeinkben minden esetben megjelenik a β paraméter, amit a megfigyelt adatainkból kell becsüljünk. Ezt a becslést tekintjük át ebben az alfejezetben. A homogén Poisson-folyamat intenzitásáról már esett szó. Szemléletesen, az intenzitás a véletlen pontok számát adja meg a vizsgált régió d-dimenziós egységtérfogatában. Legyen a B ⊂ Rd régióban az események száma N (B) (ahogy már fentebb, a Poissonfolyamatoknál definiáltuk). Ha az N (..) stacionárius akkor EN (B) = βλd (B). Így β torzítatlan becslése : N (B) , βˆ = λd (B) ahol λd a d dimenziós térfogatot jelöli [9].
2.3.2. Három Monte Carlo teszt Ahhoz, hogy elfogadjuk a CSR hipotézisünket egy X térbeli pontmintára, amely n pontból áll egy adott B ⊂ Rd (ahol d = 2, vagy 3) régióban, teszteket kell elvégezzünk. Legyen a nullhipotézisünk az, hogy a kérdéses pontminta CSR, azaz H0 = CSR. Ahhoz hogy a nullhipotézist elfogadjuk, vagy elutasítsuk a V teszt statisztika felhasználásával, α szignifikancia szinten, a következő Monte Carlo módszereket használhatjuk. Pontonkénti Monte Carlo Legyen az adott régióban a megfigyelt adatok teszt statisztikájának tapasztalati eloszlás függvénye Vˆ (r). Generáljunk s − 1 független CSR-szimulációt (ezt úgy tesszük meg, hogy megbecsüljük az X pontminta intenzitását, és ezzel az intenzitással generálunk s − 1 darab pontmintát amelyek mind a homogén Poisson-folyamat realizációi) a vizsgálandó B régióban, 14
és számoljuk ki mindegyikre a Vˆ (j) (r) (j = 2,3, . . . , s − 1) függvényt. Határozzuk meg a szimulált adatok Vˆ (j) (r) függvényei közül egy előre rögzített r értékre a minimumot és a maximumot, azaz L(r) = min Vˆ (j) (r) j
U (r) = max Vˆ (j) (r) j
Annak a valószínűsége, hogy a Vˆ (r) teszt statisztika az [L(r), U (r)] sávon kívül esik 2s , hiszen P(Vˆ (r) < L(r)) = 1 és P(Vˆ (r) > U (r)) = 1 . Így a nullhipotézist akkor utasítjuk el α = 2 s
s
s
szignifikancia szinten a fix r érték mellett, ha Vˆ (r) kivűl esik az [L(r), U (r)] sávon. Ez a módszer tehát s − 1 = 39 szimulációra α =
2 40
= 0.05 szingnifikancia szintet biztosít [21].
Szimultán Monte Carlo A pontonkénti módszernél használt, előre lerögzített r érték befolyásolhatja az [L(r), U (r)] sávot, és így a teszt kimenetelét is. A következőkben egy olyan módszert tárgyalunk, amelyhez nem kell előre lerögzített r értéket használni. Ahogy az előző módszernél, most is számítsuk Vˆ (r)-t és az s−1 szimulált CSR-hez tartozó Vˆ (j) (r) (j = 2,3, . . . , s − 1) függvényeket. Ezután számítsuk ki a teszt statisztika CSR melletti elméleti eloszlásfüggvényének és a fenti függvényeknek a maximum eltérését, azaz D(j) = max |Vˆ (j) (r) − Vpoi (r)| r
értéket, ahol Vpoi jelentse a teszt statisztika elméleti eloszlásfüggvényét, ahol j = 2, . . . , s. Minden M = s − 1 szimulált adathalmazra számítsuk ki a D(j) értéket és vegyük közülük a maximumot, amit nevezzünk Dmax -nak. Ekkor a következő felső és alsó határt kapjuk: L(r) = Vpoi (r) − Dmax U (r) = Vpoi (r) + Dmax A Vˆ (r) érték csak akkor esik kívül az [L(r), U (r)] sávon, ha a D(1) = max |Vˆ (r) − Vpoi (r)| érték nagyobb mint a Dmax érték. Ez a H0 mellett
1 M +1
r
valószínűséggel következik be, tehát
5%-os szignifikancia szinthez M = s − 1 = 19 szimuláció szükséges [21]. Rendezett Monte Carlo Alkalmazzuk a V teszt statisztikát a pontmintánkra, a kapott értéket pedig nevezzük el v1 -nek. A vizsgált B régióban generáljunk s − 1 független CSR szimulációt a homogén Poisson-folyamatból, és mindegyik szimulációra számítsuk ki a V tesztstatisztika értékeit. A kiszámított értékeket nevezzük el a következőképpen: vi : i = 2 . . . s. Jelöljük vj -vel a j. legnagyobb értéket a vi : i = 1, . . . , s értékek között. Rendezzük a kapott értékeket nagyság szerint, ekkor a következőt kapjuk: v(1) < v(2) < v(3) < . . . < v(s) . 15
Vizsgáljuk meg az alábbi valószínűséget : P (v1 = vj ) = azaz
1 s
(s − 1)! 1 = s! s
annak a valószínűsége, hogy a v1 érték a j. legnagyobb az v1 , . . . , vs között. Tehát
annak a valószínűsége, hogy v1 a rendezett v(i) értékek között a k., vagy annál nagyobb helyen áll, éppen
k s
[10], [21].
Lényegében az v1 értéket szeretnénk valamilyen értelemben elhelyezni a többi érték között úgy, hogy ha v1 a többi értékhez képest nagy (erre utal v1 helye a rendezett értékek között), akkor a CSR hipotézisünket elvetjük. Természetesen, feladattól függően akkor is elvethetjük a hipotézisünket, ha az v1 érték kicsi. A fenti módszer az első esetben adott α =
k s
szignifikancia
szinten egy jobboldalú, a második esetben baloldalú tesztet biztosít a CSR hipotézisre. Két oldalú teszthez α =
k 2s
szignifikancia szintet használunk. J. Diggle hivatkozik [10] Hope több
példájára, ahol látható, hogy s nem kell túl nagy legyen.
2.3.3. CSR tesztek A következőkben átnézzük a távolságon alapuló CSR teszteket. A továbbiakban X jelöli azokat a pontfolyamatokat, amelyek pontminták is egyben. Minden teszt esetében feltesszük, hogy a vizsgálni kívánt X pontfolyamat stacionárius. A stacionárius pontfolyamat definíciója a következő: 9. Definíció. Egy Rd -beli X pontfolyamat akkor stacionárius, ha X + v pontfolyamatnak ugyanaz az eloszlása, mint X-nek, minden v ∈ Rd -re. Legközelebbi szomszéd távolság A legközelebbi szomszéd távolság (nearest neighbour distances) [10], [21], [9] teszt-statisztika n eseményre, egy B ⊂ Rd régióban, az X stacionárius pontfolyamat xi pontjáhozhoz legközelebbi pont távolságát méri, azaz a ti = min ||xi −xj ||, vagyis az összes xj ∈ X, (i 6= j) pontpár j6=i
távolságának a minimumát. Ennek a teszt statisztikának az eloszlásfüggvénye a következő: 10. Állítás. A legközelebbi szomszéd távolság eloszlásfüggvénye: (2.1)
G(r) = P (dmin (x, u) < r|x, u ∈ X)
ahol x egy tetszőleges esemény, és dmin (x, u) az x eseményponttól a hozzá legközelebbi esemény távolságát jelöli az X pontmintában. A fenti G(r) függvénynek a tapasztalati eloszlásfüggvénye a következő: P ˆ 1 (r) = i 1{ti ≤ r} (2.2) G n
16
A G(r) elméleti eloszlásfüggvénye két dimenzióban G2poi (r) = 1 − exp(−βπr2 ),
(2.3) három dimenzióban pedig: (2.4)
4 G3poi (r) = 1 − exp(− βπr3 ), 3
ahol a β intenzitást a fenti becsléssel kapjuk meg. A két függvényt mindhárom Monte Carlo módszerrel tesztelhetjük. A szimultán és a ponˆ 1 (r) > Gpoi (r) akkor a Poisson folyamatnál rövidebb utak vannak a tonkénti MC-re, ha G ˆ 1 (r) < Gpoi (r) pedig szabályos pontminpontfolyamatban, ami klaszterezett pontmintára a G tára utal. Rendezett MC-re a következő értékeket definiálhatjuk : 1. ui = t¯, a ti -k középértékét, R 2 ˆ2 ˆ (r) − G ¯ 2 (r))2 dr, ahol G ¯ 2 (r) = (s − 1)−1 P G 2. ui = (G i i i i6=j i (r) Pont és a hozzá legközelebbi esemény távolsága A következő statisztikánk a szakirodalomban a "pont és a hozzá legközelebbi esemény távolsága" (point to nearest event distances), vagy az "üres tér" távolságok (empty space distances, és inenntől F függvényként hivatkozunk rá) néven jelenik meg [10], [21], [9]. Ez a módszer a következőképpen működik. Vegyünk fel m darab, véletlenszerű pontot a megfigyelt B ⊂ Rd régióban (ezt véletlen pontgenerálással érhetjük el), és jelölje di az m elemszámú pontmintából tetszőlegesen választott vi referencia pont és a hozzá legközelebb lévő pontfolyamatbeli x pontnak a távolságát. Másképp, legyenek X stacionárius pontfolyamat pontjai az xj pontok és legyen vi egy tetszőleges referencia pont. Ekkor a pont és a hozzá legközelebb eső esemény távolsága formálisan a következőt jelenti : dmin (vi , xj ) = min(||vi − xj || : xj ∈ X). i
11. Állítás. Az F eloszlásfüggvénye: (2.5)
F (r) = P (d(v, X(x)) ≤ r),
ahol v egy tetszőleges referencia pont. Az Fˆ tapasztalati eloszlásfüggvény egy v1 , . . . , vm referencia pontokból álló B régióban: P 1{dmin (vi , xj ) ≤ r} ˆ . (2.6) F (r) = i m Az F (r) elméleti eloszlásfüggvény két dimenzióban (2.7)
2 Fpoi (r) = 1 − exp(−βπr2 ),
17
három dimenzióban pedig: (2.8)
4 3 Fpoi (r) = 1 − exp(− βπr3 ). 3
CSR hipotézisünk eldöntéséhez itt is használhatjuk mind a három MC-módszert. Mint fentebb a G függvénynél, itt is érvényes az, ha Fˆ (r) > Fpoi (r) akkor feltételezhetően szabályos pontmintával van dolgunk, és ha Fˆ (r) < Fpoi (r), akkor pedig klaszterezettel. A rendezett Monte Carlohoz használhatjuk az R – ui = (Fˆi (r) − F¯i (r))2 dr ¯ i (r) függvényt. értéket, ahol F¯i (r) függvényt ugyanúgy kapjuk minta G Ripley féle K függvény A K függvény az X pontfolyamat egy xi pontjától a legfeljebb r távolságban lévő pontok várható számának és a β intenzitásnak a hányadosát adja meg [10], [21], [9]. Tegyük fel, hogy az X pontfolyamatra a szokásos tulajdonságunk a stacionaritás mellett, az izotrópia (elforgatás stacionaritás) is teljesül. Ekkor a K függvény:
(2.9)
K(r) =
E(Nx (r)) , β
ahol Nx (r) az X pontfolyamat egy tetszőleges pontjától r távolságban lévő pontok számát adja meg. Ez, a fenti feltételek mellett, homogén Poisson-folyamatra két dimenzióban: (2.10)
2 Kpoi (r) = πr2 ,
háromdimenzióban pedig a: (2.11)
4 3 Kpoi (r) = πr3 3
elméleti eloszlásfüggvényeket adja. ˆ többféleképpen is megkaphatjuk, de itt A K függvény tapasztalati eloszlásfüggvényét, K-t csak a következő módszert vizsgáljuk. Legyen n eseményünk egy B régióban, és vegyük az X pontfolyamat xi és xj (ahol i 6= j) pontjai között a tij = ||xi − xj || távolságot. Ekkor a tapasztalati eloszlásfüggvény: P P x ∈B i6=j 1{||xi − xj || ≤ r} i ˆ = (2.12) K . ˆ d (B) βλ ˆ ≥ Kpoi , akkor Itt is használhatjuk a szimultán, vagy a pontonkénti Monte Carlot. Ha K ˆ ≤ Kpoi , akkor szabályos pontmintát feltételezhetünk. klaszterezett, ha K
18
A χ2 módszer Ezt a módszert az előzőekben tárgyalt távolságon alapuló módszerek alternatívájaként mutatjuk be. Legyen adva egy B kockában egy n elemű X pontminta. Osszuk fel B-t m darab C1 , . . . Cm egyenlő nagyságú, kisebb kockára, ahol m egy általunk választott szám. Legyen i1 , . . . , im a C1 , . . . Cm kis kockákba eső események száma, azaz ij = #(Cj ∩ X). Jelölje továbbá n0 , n1 , . . . , nm azoknak a C1 , . . . , Cm kis kockáknak a számát, amelyekbe rendre 0,1, . . . , m esemény jutott, azaz n0 = #(ij = 0), n1 = #(ij = 1), . . . , nm = #(ij = m). P n Becsüljük a Poisson-folyamat β várhatóértékét : βˆ = m , ahol i ni = n. Ekkor meghatározhatjuk a pˆ0 , . . . , pˆm értékeket, ahol pˆk =
ˆ βˆk −λ . k! e
m X (nl − nˆ pl )2 l=0
nˆ pl
Így a χ2 -statisztika a következő lesz :
∼ χ2m−1
Ez a tesztelési módszer függ a felosztás finomságától, így több felosztásra is tesztelni fogjuk P ˆ2 2 (ni −β) a hipotézisünket. Említsük meg azt is, hogy a σˆ arány, ahol σ 2 = m i=1 m−1 , Poissonβ
eloszlás esetén egyenlő eggyel, így ezt tulajdonságot is vizsgálhatjuk az elemzéseinknél [21].
19
2.4. Megjegyzések a fejezethez Ebben az alfejezetben a fentebb elhangzottakhoz fűzünk hozzá néhány állítást, meggondolást, érdekességet.
2.4.1. F, G, K függvények Az F és a G függvények elméleti eloszlása Az alábbiakban bemutatjuk a G elméleti eloszlás függvényének levezetését, a kalkulus eszközeit használva. Ahogy az előző fejezetben láttuk: G(r) = P (dmin (x, u) < r|x, u ∈ X). Tegyük fel, hogy a vizsgálandó X pontfolyamat n darab pontja a sík egy B tartományában helyezkedik el. Legyen x egy tetszőleges pontja az X pontfolyamatunknak. A többdimenziós Poisson folyamat egyik, már az előző fejezetekben tárgyalt tulajdonsága : ha N (B) = n akkor az n darab esemény független minta a többdimenziós egyenletes eloszlásból. Ekkor annak a valószínűsége, hogy az x pont egy r sugarú környezetében található egy x-től különböző pontfolyamatbeli pont: (2.13)
r2 π . λ2 (B)
A fenti összefüggésből könnyen megkaphatjuk a komplementer esemény valószínűségét, ami nem más mint : (2.14)
1−
r2 π , λ2 (B)
azaz ennyi annak a valószínűsége, hogy az x esemény r sugarú környezetében nem található x-től különböző esemény. Tudjuk, hogy a vizsgált n eseményünk független, ezért az összes pontra meghatározva az utóbbi valószínűséget, a következőt kapjuk: r2 π n (2.15) 1− . λ2 (B) Korábban azt is láttuk, hogy a β intenzitásra a
N (B) λ2 (B)
becslést használtuk. Ebben az esetben
N (B) = n, a feltétel miatt. Ezt észben tartva szorozzuk be a (2.15) belső függvényét β-val, majd vizsgáljuk meg az így keletkezett összetett függvény határértékét: r2 πβ n 2 (2.16) lim 1 − = e−βr π n→∞ n
20
n λ2 (B)
=
A fenti (2.16) összefüggésből arra következtetünk, hogy elég nagy n-re (2.15)-nek jó közelítése a (2.16), azaz (2.17)
r2 π n 2 1− ∼ e−βr π . λ2 (B)
Tehát elég nagy n esetén annak a valószínűsége, hogy az X pontfolyamat pontjainak az r környezetében található pont: (2.18)
2
1 − e−βr π ,
ami éppen a G és F függvények elméleti eloszlása két dimenzióban. A fenti gondolatmenet pontosan ugyanez magasabb dimenziókban is. A K függvény eloszlása Láttuk, hogy a Ripley-féle K függvény eloszlása: (2.19)
K(r) =
ENx (r) β
alakú, ahol a számláló egy pontfolyamatbeli x esemény r sugarú környezetében a várható események számát, a nevező pedig a homogén Poisson-folyamat várható érték jelenti a vizsgált régióban. Két dimenzióban a számláló a βr2 π alakot veszi fel, és így a hányadost β-val egyszerűsítve megkapjuk a keresett függvényt. (2.20)
K(r) =
βr2 π = r2 π β
Ugyanez a gondolatmenet alkalmazható magasabb dimenziókban is.
2.4.2. Határ problémák Térbeli pontminták elemzése során akkor merülnek fel a határ problémák, amikor a vizsgálandó régiónk részhalmaza egy olyan nagyobb régiónak, ahol nincsen adatunk a pontfolyamat többi pontjáról. A probléma abból fakad, hogy nehéz megállapítani, hogy a vizsgált régión kívül eső események milyen kapcsolatban állnak a vizsgált régióban lévő eseményekkel [21]. Ezt 2.5 ábra is szemlélteti, ahol a sárgával megjelölt tartomány a vizsgált régió, a régióban lévő események pedig feketével vannak színezve. Példaként vegyük azt az esetet, amikor a legközelebbi szomszéd statisztikát használjuk. Ekkor nem tudjuk eldönteni biztosan, hogy egy pontfolyamatbeli esemény egy a vizsgált régióban lévő eseményhez, vagy egy a vizsgált régión kívüli eseményhez van-e közelebb ? Az ilyen problémák megoldására határ korrekciót használunk. A dolgozatban a vizsgált távolságon alapuló statisztikák mindegyikére határ korrekciót kell alkalmazzunk minden olyan esetben, ahol ezt a feladat körülményei megkövetelik. 21
2.5. ábra. Határ probléma [21]
2.4.3. Legközelebbipont algoritmusok A legközelebbi pontok megtalálásának a feladata, a számítógépes geometriának egy intenzíven kutatott területe. Ebben az alfejezetben a teljesség igénye nélkül, említés szintjén bemutatunk néhány eljárást a legközelebbi pontpárok megtalálására. Nyers erő Tegyük fel, hogy adott a d-dimenziós térben egy V ponthalmaz, amelyben n pont található. Ezen a téren legyen adott egy Lq metrika. A legközelebbi pont probléma a következőképpen írható le: minden p ∈ V pontra találjuk meg V − {p} halmazban a p-hez legközelebbi szomszédot, az Lq metrikát felhasználva. Az egyszerűség kedvéért vizsgáljuk azt az esetet, amikor a ponthalmazunk a d-dimenziós 1 P 2 2, euklideszi térben található, és Lq = L2 = i |xi − yi | A legegyszerűbb módszer a feladat megoldására az úgynevezett nyers erő módszer, amely egy tetszőleges pi pontra kiszámítja az összes p1 , p2 , . . . , pi−1 , pi+1 , . . . , pn pont euklideszi távolságát, majd a kiszámított értékek közül, kiválasztja a minimumot. Ezzel a módszerrel az összes legközelebbi pontpár megtalálásához O(n2 ) futási idő szükséges. Voronoi-cellák Legyen a fenti feltétekkel megfogalmazott, szabálytalan elrendezésű ponthalmazunk a síkon. Vizsgáljuk meg a következő állításokat [5]: 1. Minden pi ∈ V pont köré szerkeszthető olyan sokszög, amelynek a belső pontjai közelebb esnek pi -hez, mint az összes többi V − {pi } ponthoz. 2. Az ilyen sokszögek konvexek, és folytonosan kitöltik a síkot.
22
A fenti meghatározásból következik, hogy a sokszög oldalai merőlegesek a pontokat összekötő egyenesekre és felezik azokat. Nem mindegyik felező merőleges lesz a sokszög része, csak azok, amelyek zárt metsződéseiből a konvex sokszög létre jön. A sokszög oldali egyben meghatározzák a szomszédos pontokat is, hiszen azok a pontok befolyásolják a sokszög alakulását.
2.6. ábra. Voronoi felbontás A Voronoi-cellák meghatározása után az egyes pontok legközelebbi szomszédait már egyszerű megtalálni, hiszen a szomszédok már adottak, így a futási idő lényeges részét a cellák kiszámítása viszi el. Két dimenzióban a cellák kialakítása O(n log n) időben lehetséges [6]. Magasabb dimenziókban ez a futási idő lényegesen rosszabb. Legközelebbi pontok O(n log n) időben Ez az alfejezet egy O(n log n) futásidejű eljárást ismertet, rögzített d-dimenziós térben Lq metrika mellett. Ez a módszer egy a J. L. Bentley által javasolt, szintén optimális módszernek egy módosítása. A cikk [14] írójának álláspontja és tapasztalata szerint az átlagos gyakorlati esetekben ez egy gyorsabban működő, determinisztikus algoritmus. Az algoritmus lényege a következőkben foglalható össze. Az adott d dimenziós ponthalmazt egy izotetikus kocka halmaz struktúrában kezeljük, amivel párhuzamosan minden egyes kockához vezetjük a "szomszédsági" információt leíró adatszerkezetet. A kockákat bizonyos rendszer szerint kisebb kockákra osztjuk, amivel egyidejűleg szétosztjuk az eredeti kockában lévő próbapontokat tartalmazás szerint az új, kisebb kockákba, és aktualizáljuk a szomszédsági kapcsolatokat. Az eljárást addig ismételjük, amíg minden egyes megmaradó kocka egyetlen próbaponttá "zsugorodik". Ekkor az eljárás terminál, és minden eredménykockához (próbaponthoz) a szomszédsági struktúrából közvetlenül leolvasható a hozzá legközelebbi pontok halmaza. Az algoritmus komplett és precíz leírása nagyon terjedelmes lenne. Számos algoritmuselméleti trükk alkalmazásával és diszkrét geometriai lemma bizonyításával igazolja, hogy a módszer futásideje O(n log n), és tárhely igénye pedig O(n). A cikk továbbá igazolja, hogy ennél jobb futási időt elméletileg sem lehet elérni az általános problémára.
23
3. fejezet
Elemzések Ebben a fejezetben az első részben tárgyalt, csillagkeletkezésben aktív galaxisok eloszlását fogjuk elemezni a második fejezetben tárgyalt néhány módszerrel.
3.1. A vizsgálni kívánt adatok leírása A vizsgált adatmátrixunk 10000 sorból, és 4 oszlopból áll. Mindegyik sor egy galaxisnak felel meg, az oszlopok pedig a galaxisok következő paramétereit tárolják: – galaxyID : a galaxis nevét, – (x,y,z) : a galaxis koordinátáit. A galxisok ∼ 500 × 500 × 500 térfogatú kockában találhatóak, ez lesz a vizsgálandó régiónk. Az adatainkat a Millennium szimuláció adatbázisából, egyszerű SQL parancsokkal kérdeztük le és mentettük el egy .csv kiterjesztésű fájlba. Az adatok R programba való beolvasáshoz a read.csv() parancsot használtuk, majd a spatstat package betöltése után a galaxisok nevű, pp3 háromdimenziós pontminták tárolására használt objektumba mentettük el.
3.2. Távolságon alapuló tesztek 3.2.1. A G függvény Az R programcsomagban a legközelebbi szomszéd teszt statisztikát három dimenzióban a spatstat : :G3est() paranccsal érhetjük el. Ha a G statisztikát a szimultán Monte Carloval szeretnék használni, akkor a következő kódot gépeljük be : – V <- envelope(galaxisok,G3est,global=TRUE,nsim=39,nrank=1),
24
ahol a global=TRUE paraméterrel a szimultán Monte Carloval való számolást, az nsim=39 argumentummal az adott régióban létrehozott 39 CSR szimulációt, az nrank=1 paraméterrel pedig a maximális és minimális sávokkal való számolást állítottuk be. A szignifikancia szint a szimultán MC fejezetben tárgyaltak alapján: α =
1 40 .
A CSR
hipotézisünket akkor vetjük el α szignifikancia szinten, ha a tapasztalati G függvény bármilyen r értékre kilép a szimultán MC-vel meghatározott sávokon. Ezt vizsgálhatjuk meg a 3.1 ábrán.
3.1. ábra. A G függvény, szimultán MC-vel Bár a nagymennyiségű megfigyelt esemény ellenére is lehet látni, hogy a G3obs -nak nevezett tapasztalati függvény kilép a megadott G3hi sávon, ezt, a biztonság kedvéért, a következő programmal is leellenőrizhetjük. logikaiVektorG3 <- G3$hi>G3$obs table(logikaiVektorG3)["FALSE"]>0 Itt az első sorban létrehoztunk egy olyan logikai vektort, amely TRUE értékeket fog tárolni azokon a helyeken ahol a G függvény nem lépett ki, és FALSE értékeket azokon a helyeken ahol pedig kilépett a megadott sávból. A második sorban meghatározzuk, hogy milyen gyakran jelent meg FALSE érték ebben a logikai vektorban, és ha a kapott szám nagyobb mint nulla, akkor a hipotézisünket, α szignifikancia szinten elvetjük. A programot lefuttatva nullánál több FALSE értéket kaptunk vissza, így a CSR hipotézisünket elvetjük.
3.2.2. Az F függvény Hasonlóan a G függvényhez, az "üres tér" teszt statisztikát a spatstat : :F3est() paranccsal érhetjük el. 39 CSR szimulációval, szimultán Monte Carloval meghatározva a sávokat, a kö25
vetkező ábrát kapjuk :
3.2. ábra. Az F függvény, szimultán MC-vel 3 -nak nevezett tapasztalati függvény kilép az F 3 Ennél az ábránál is látszik, hogy az Fobs lo
sávon. Ezt újra ellenőrizhetjük a G függvénynél bevált logikai vektort használó programunkkal. Kimenetnek egy nullától különböző számot kaptunk, így a CSR hipotézisünket ismét el kell vessük α szignifikancia szinten.
3.2.3. A K függvény A fenti statisztikákhoz hasonlóan, a K függvényt a spatstat : :F3est() paranccsal hívhatjuk meg. Szimultán MC-t használva, 39 CSR szimulációval a 3.3 ábrát kapjuk. 3 és K 3 sávokon. A fenti programot 3 függvény nem lép ki a Klo Látható, hogy a Kobs hi
lefuttatva, a kapott kimenet alátámasztja ezt az észrevételt. A kapott eredmények fényében, a CSR hipotézisünket nem vetjük el.
26
3.3. ábra. A K függvény, szimultán MC-vel
3.3. A χ2 -próba A χ2 -statisztikát több felosztásra is el fogjuk végezni. A kapott eredményeket az alábbi táblázat tartalmazza: Felosztás finomsága
α
df
χ2krit
χ2kapott
Ítélet
106
0,05
3
7,8174
283,43
elutasítjuk H0 -t
0,05
4
9,4877
316,21
elutasítjuk H0 -t
0,05
8
15,5073
956.84
elutasítjuk H0 -t
0,05
11
19,6751
2623,7
elutasítjuk H0 -t
1,25 · 105 1,5625 ·
104
8 · 103
A fenti táblázatban a H0 = (a pontmintánk CSR) hipotézist vizsgáltuk több finomodó felosztás mellett. Mindegyik felosztásra elvettük a H0 hipotézist, így a χ2 -próbát használva, kimutattuk, hogy a pontmintánk nem CSR.
27
4. fejezet
Összefoglalás és további kérdések A feladatomnak csillagkeletkezésben aktív galaxisok eloszlásának a vizsgálatát választottam. Ez a probléma érdekes kérdés, tudva, hogy a gammakitörések eloszlása nem véletlenszerű, és hogy a GRB-k csillagkeletkezésben aktív galaxisokban jönnek létre. A valószínűségszámítás és a statisztika eszköztárát felhasználva, az elemzések és vizsgálatok végén arra a konklúzióra jutottam, hogy a csillagkeletkezésben aktív galaxisok eloszlása nem véletlenszerű. Nyitott kérdés maradt azonban, hogy ha az eloszlás nem véletlenszerű, akkor vajon milyen szerkezetek fedezhetőek fel, és milyen matematikai modellek ilesszkednek a legjobban a megtalált struktúrákra? Ezen felül nyitott kérdés maradt még az is, hogy ha a galaxisok koordinátáin kívül más paramétereket is figyelembe veszünk, mint például a galaxis tömegét, a galaxisban megtalálható feketelyuk tömegét, a sötétanyag sűrűségét a galaxis körül, stb. . ., akkor vajon ugyanerre az eredményre jutottunk volna-e ?
28
Irodalomjegyzék [1] http://photojournal.jpl.nasa.gov/jpeg/PIA16876.jpg. v, 3 [2] http://www.csillagaszat.hu/wp-content/uploads/2014/02/20140201_ kvazarhalmaz_fig1.png. v, 4 [3] http://www.mpa-garching.mpg.de/galform/virgo/millennium/seqB_063a_half. jpg. v, 5 [4] http://www.mpa-garching.mpg.de/galform/virgo/millennium/seqF_063a_half. jpg. v, 5 [5] http://www.agt.bme.hu/tutor_h/terinfor/t27.htm. 22 [6] http://hu.wikipedia.org/wiki/Voronoj-cella. 23 [7] G. Ákos. A xx. századi csillagászat csonka koronája : gamma-kitörések kutatása. 3, 4 [8] A. Baddeley, I. Bárány, and R. Schneider. Spatial point processes and their applications. Stochastic Geometry : Lectures given at the CIME Summer School held in Martina Franca, Italy, September 13–18, 2004, pages 1–75, 2007. 9, 10, 11, 12, 13 [9] A. Baddeley, R. Moyeed, C. Howard, and A. Boyde. Analysis of a three-dimensional point pattern with replication. Applied Statistics, pages 641–668, 1993. 14, 16, 17, 18 [10] P. J. Diggle et al. Statistical analysis of spatial point patterns. Academic press, 1983. v, 11, 13, 16, 17, 18 [11] I. Horváth. Gammakitörések mutatják a világegyetem legnagyobb szerkezetét. http:// www.csillagaszat.hu/hirek/ko-korai-vilagegyetem/ko-univerzum-szerkezete/ gamma-kitoresek-mutatjak-a-vilagegyetem-legnagyobb-szerkezetet/. 4 [12] I. Horváth, J. Hakkila, and Z. Bagoly. Possible structure in the grb sky distribution at redshift two. Astronomy & Astrophysics, 561 :L12, 2014. 1, 5 [13] V. Springel, S. D. White, A. Jenkins, C. S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly, et al. Simulating the joint evolution of quasars, galaxies and their large-scale distribution. arXiv preprint astro-ph/0504097, 2005. 5 [14] P. M. Vaidya. Ano (n logn) algorithm for the all-nearest-neighbors problem. Discrete & Computational Geometry, 4(1) :101–115, 1989. 23 [15] Wikipedia. Big bang. http://en.wikipedia.org/wiki/Big_Bang. 2 [16] Wikipedia. Cosmological principle. http://en.wikipedia.org/wiki/Cosmological_ principle. 1 29
[17] Wikipedia. Gamma ray. http://en.wikipedia.org/wiki/Gamma_ray. 3 [18] Wikipedia. Gamma-ray burst. http://en.wikipedia.org/wiki/Gamma-ray_burst. 1 [19] Wikipedia. Millennium run. http://en.wikipedia.org/wiki/Millennium_Run. 5 [20] Wikipedia. Ösrobbanás. 2 [21] V. Ying. Methods for spatial analysis on a network. Master’s thesis, UCLA. v, 15, 16, 17, 18, 19, 21, 22 A megjelölt weboldalakat, 2015 máj 25-én leellenőriztem.
30
NYILATKOZAT
Név: Daróczy-Kiss Endre ELTE Természettudományi Kar, szak: Elemző matematikus Neptun azonosító: D71ZRG Szakdolgozat címe: Galaxis-eloszlások vizsgálata
A szakdolgozat szerzőjeként fegyelmi felelősségem tudatában kijelentem, hogy a dolgozatom önálló munkám eredménye, saját szellemi termékem, abban a hivatkozások és idézések standard szabályait következetesen alkalmaztam, mások által írt részeket a megfelelő idézés nélkül nem használtam fel.
Budapest, 2015.05.30
_______________________________ hallgató aláírása
