preizkušanje najbližji sosed box-mueller neodvisnost slu indikatorji monte č semivariogram ajno domnev križni variogram polje carlo gumbel fakulteta za gradbeništvo in geodezijo navadno krigiranje rezidual univerza v ljubljani weibull PDF poisson 2023 frechet gostota prag preizkus hipoteza moran robustna prostorska geostatistika kovarinca igla prejem s trendom ekstremgeary odstopanje buffonova zavrnitev krigiranje inverznavezanikrigiranje CDF kovariogram porazdelitena osamelec metoda razsevni krigevarianca statistika graf funkcija vzorec jura najmanjši kvadrati goran turk regresija kvadrati variogram cholesky GORAN TURK Prostorska statistika UNIVERZITETNI U ČBENIK LJUBLJANA, OKTOBER 2023 PROSTORSKA STATISTIKA Univerzitetni učbenik Avtor: Goran Turk Recenzenta: Gregor Dolinar, Samo Drobne Jezikovni pregled: Nina Ambrožič Naslovnica: Goran Turk Oblikovanje in prelom: Goran Turk Založnik: Založba Univerze v Ljubljani (University of Ljubljana Press) Za založbo: Gregor Majdič, rektor Univerze v Ljubljani Izdajatelj: Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo Za izdajatelja: Violeta Bokan Bosiljkov, dekanja Fakultete za gradbeništvo in geodezijo Ljubljana, 2023 Prva e-izdaja. Publikacija je v digitalni obliki prosto dostopna na: https://ebooks.uni-lj.si DOI: 10.15292/9789612971991 Kataložni zapis o publikaicji (CIP) pripravili v Narodni in univerzitetni knjižnici v Ljubljani COBISS.SI-ID 171155459 ISBN 978-961-297-199-1 (PDF) To delo je ponujeno pod licenco Creative Commons Priznanje avtorstva-Nekomercialno-Brez predelav 4.0 Mednarodna licenca. This work is licenced under a Creative Commons Attribution- NonCommercial-NoDerivs 4.0 International licence. Predgovor Prostorska statistika je razširitev običajne statistike, pri kateri nimamo podatka o legi ali času izvedbe poskusa ali pa nas to ne zanima. Pri prostorski statistiki obravnavamo slučajne spremenljivke z dodatno informacijo o legi ali času, za kateri smo dobili posamezen element vzorca. Takim slučajnim spremenljivkam pravimo slučajna polja, če imamo podatke o legi v prostoru, ali slučajni procesi, če poznamo podatke o času ali legi v prostoru in času. Učbenik vključuje kratko uvodno poglavje, ki opredeli osnovne značilnosti prostorske statistike, prikaže razvoj oziroma zgodovino te relativno nove veje uporabne statistike in pojasni, s katerimi vrstami spremenljivk in podatkov se prostorska statistika ukvarja. V tem poglavju je tudi seznam različnih geostatističnih podatkov, ki jih študenti lahko uporabijo za vaje in preizkušanje različnih metod prostorske statistike. Večino teh podatkov sem pridobil na svetovnem spletu, primer podatkov iz Slovenije pa sta prispevala raziskovalca Kmetijskega inštituta Slovenije, Borut Vrščaj in Janez Bergant, za kar se jima iskreno zahvaljujem. Za razumevanje prostorske statistike je nujno osnovno predznanje verjetnostnega računa in statistike. To predvsem vključuje poznavanje in razumevanje osnovnih pojmov, kot so slučajna spremenljivka, porazdelitev slučajne spremenljivke, porazdelitvena funkcija, slučajni vektor, momenti slučajnih spremenljivk in vektorjev (pričakovana vrednost, varianca, kovarianca), funkcije slučajnih spremenljivk in momenti funkcij slučajnih spremenljivk, osnovne statistike (povprečje, varianca vzorca, kovarianca vzorca), preizkušanje domnev ter regresija. Zato je v učbeniku kratko drugo poglavje, kjer so na kratko našteti in definirani ti osnovni pojmi. V zaključku tega poglavja je definirana tudi skalarna in vektorska slučajna funkcija. V poglavju Geostatistika so opisane osnovne metode za analizo slučajnih polj in procesov. To vključuje definicijo razsevnih grafov, kovariančnih funkcij, variograma in dveh pogosto uporabljenih preizkusov domnev o lastnostih slučajnih polj (Moranov indeks in Gearyjevo razmerje). Krigiranje je pomembno regresijsko orodje, zato ga tu opišemo v samostojnem poglavju. Opišemo tri osnovne metode krigiranja: preprosto, običajno in krigiranje s trendom. Kot alternativa krigiranju lahko uporabimo tudi prostorsko regresijo, ki je pravzaprav običajna linearna regresija z utežmi, ki so določene tako, da bližnje točke močneje vplivajo na vrednost regresije. Postopki določitve uteži in parametrov prostorske linearne regresije so opisani v poglavju Prostorsko utežena re- gresija. V poglavju Prostorski vzorci analiziramo razpored točk v prostoru, same vrednosti polja pa nas tu ne zanimajo. Opišemo dve skupini metod: metoda kvadratov in metoda najbližjega soseda. Zadnje poglavje na videz ne sodi v učbenik Prostorska statistika, saj je metoda Monte Carlo in generiranje vzorcev slučajnih spremenljivk bolj splošno področje. Kljub temu pa predstavlja osnovni gradnik pri generiranju slučajnih polj in procesov in zato kot zadnje poglavje ustrezno zaključuje učbenik. Pisanje učbenika je dolgotrajen postopek, ki se je pri meni vedno začel s pisanjem zapiskov pred in po predavanjih. Po nekaj letih je nastal osnutek učbenika, ki sem ga nato uredil in pripravil v obliki, ki je bila že precej podobna zdajšnji verziji. V nekaj naslednjih letih smo jaz in moji študenti učbenik upo-rabljali in ob tem odkrivali različne napake, pomanjkljivosti in možnosti za izboljšanje. Študentom, ki so s pripombami aktivno sodelovali pri tem, se na tem mestu lepo zahvaljujem. V zaključni fazi priprave učbenika sta učbenik skrbno prebrala recenzenta Gregor Dolinar in Samo Drobne. Za njune pripombe, ki nedvomno izboljšujejo učbenik, se jima iskreno zahvaljujem. Posebej pa se moram zahvaliti še Nini Ambrožič, ki je učbenik zelo natančno lektorirala. Že na začetku pisanja sem se odločil, da učbenik ne bo šel v klasični tisk, temveč bo na voljo prvenstveno v digitalni obliki. Zato učbenik vključuje povezave na različne spletne strani z dodatnimi informacijami, vključene pa so tudi kratke animacije, ki še dodatno pojasnijo vplive različnih parametrov. Učbenik je v prvi vrsti namenjen študentom druge stopnje pri pouku Prostorske statistike, lahko pa tudi kot osnovno gradivo študentom tretje stopnje, ki jih zanima oziroma želijo uporabiti metode te veje statistike. Upam, da vam bo učbenik pomagal do znanja, ki si ga želite pridobiti. Goran Turk Ljubljana, jesen 2023 Vsebina Predgovor i Vsebina iii 1 Uvod 1 1.1 Zgodovina prostorske statistike . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Podatki za prostorsko statistiko . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Vrste obravnavnih slučajnih spremenljivk . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3.1 Geostatistični podatki . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.2 Podatki v mreži . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.3 Prostorski točkovni vzorci . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Slučajne funkcije, polja in procesi 9 2.1 Slučajne spremenljivke . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Slučajni vektorji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Momenti slučajnih spremenljivk in vektorjev . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Slučajne funkcije . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 Vektorske slučajne funkcije . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3 Geostatistika 23 3.1 Geostatistični podatki . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 Priprava podatkov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Razsevni graf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 Kovariančna funkcija . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 Variogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.6 Določitev empirične funkcije variograma . . . . . . . . . . . . . . . . . . . . . . . . . 31 iv Vsebina 3.7 Geostatistična analiza z indikatorji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.8 Preizkušanje domnev o prostorski povezanosti . . . . . . . . . . . . . . . . . . . . . . . 35 3.8.1 Moranov indeks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.8.2 Gearyjevo razmerje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4 Krigiranje 45 4.1 Preprosto krigiranje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2 Običajno krigiranje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.3 Krigiranje s trendom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 Primeri krigiranja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.4.1 Preprosto krigiranje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.4.2 Običajno krigiranje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.4.3 Krigiranje s trendom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5 Prostorsko utežena regresija 61 5.1 Linearna regresija . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2 Utežena linearna regresija . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2.1 Določitev uteži . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2.2 Izbira funkcije za določitev uteži in kritične razdalje . . . . . . . . . . . . . . . 65 5.2.3 Primeri . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 6 Prostorski vzorci 69 6.1 Metoda kvadratov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2 Primeri – metoda kvadratov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.3 Pomanjkljivosti metode kvadratov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.4 Metoda najbližjega soseda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.5 Primeri – metoda najbližjega soseda . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7 Metoda Monte Carlo 87 7.1 Generiranje vzorca slučajne spremenljivke . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.1.1 Inverzna metoda generiranja vzorca . . . . . . . . . . . . . . . . . . . . . . . . 89 7.1.2 Generiranje vzorca normalno porazdeljenih slučajnih spremenljivk . . . . . . . . 91 7.1.3 Metoda sprejema/zavrnitve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 7.1.4 Primeri . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.2 Generiranje vzorca slučajnega vektorja . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.2.1 Neodvisne slučajne spremenljivke – inverzna metoda . . . . . . . . . . . . . . . 98 Vsebina v 7.2.2 Odvisne slučajne spremenljivke – inverzna metoda . . . . . . . . . . . . . . . . 100 7.2.3 Korelirane normalno porazdeljene spremenljivke . . . . . . . . . . . . . . . . . 101 7.2.4 Primeri . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 8 Literatura 113 Stvarno kazalo 115 1 Uvod 1.1 Zgodovina prostorske statistike Prostorska statistika se je podobno kot druge statistične metode razvijala vzporedno s potrebami pri različnih aplikativnih problemih. Elemente prostorske statistike najdemo v 17. stoletju, ko so naredili prvo karto vetrov na Zemlji (Halley1, 1686). Znan je problem Buffonove igle, ki ga prav tako lahko sma-tramo za enega od začetkov prostorske statistike (Comte de Buffon2, 1780). Slučajni procesi (ali polja), kot osnova prostorske statistike, so se prvič uporabili v najpreprostejši obliki Poissonovega procesa pri izračunu srednje dolžine poti molekule v plinu (Clausius3, 1858). Poissonov proces še danes predstavlja osnovno merilo o enakomerni in neodvisni, a slučajni razporeditvi elementov v prostoru ali času. Že v začetku 20. stoletja so določili porazdelitev razdalje med poljubno točko v Poissonovem procesu in njeno najbližjo sosedo (Hertz4, 1909), kar je še danes osnova za nekatere metode prepoznavanja prostorskih vzorcev. Pomemben napredek v razvoju vseh statističnih metod, tudi prostorske statistike, se je odvijal v prvi polovici 20. stoletja v Angliji. Gosset5 in Fisher6 sta uporabljala metode prostorske statistike in o njih pisala, prvi v pivovarni Guinness (pod psevdenimom Student), drugi na posestvu Rothamsted Experi-mental Station v Hertforshiru. Fisher je precejšnji del svojega časa namenil statistični analizi podatkov s kmetijskih polj. Naslednji pomemben napredek se je zgodil v povezavi z geološkimi raziskavami za potrebe rudarjenja v Južni Afriki (Krige7, 1951) in v Franciji (Matheron8, 1954–64) v sredini 20. stoletja. V tem obdobju je bil pomemben tudi prispevek švedskega avtorja Mat´erna9 (1960), ki je deloval na Gozdarski fakul-teti Švedske univerze za kmetijske znanosti. Z uporabo modela zveznega prostorskega spreminjanja so postavili osnove za Gaussov10 prostorski slučajni proces, na osnovi katerega temelji pomembna veja 1 Edmond Halley, angleški astronom in matematik, znan po izračunu orbite kometa (Halleyev komet) (1656–1742). 2 Georges Louis Leclerc Comte de Buffon, francoski naravoslovec in matematik (1707–1788). 3 Rudolf Clausius, nemški matematik in fizik (1822–1888). 4 Heinrich Hertz, nemški fizik (1857–1894). 5 William Sealy Gosset, angleški statistik (1876–1937). 6 Sir Ronald Aylmer Fisher, angleški statistik in genetik (1890–1962). 7 Danie Gerhardus Krige, južnoafriški rudarski inženir, avtor postopka aproksimacije kriging (1919–2013). 8 Georges Francois Paul Marie Metheron, francoski statistik in geolog (1930–2000). 9 Bertil Matérn, švedski statistik (1907–2007). 10 Johann Friedrich Carl Gauss, nemški matematik (1777–1855). 2 1 Uvod prostorske statistike: geostatisika. Danes se prostorska statistika uporablja v različnih aplikacijah: • epidemiološke študije, • analize razvoja bolezni, • prostorsko načrtovanje, • določanje kakovosti rudnine v nahajališčih, • določanje globine in debeline geoloških plasti, • določanje prepustnosti in poroznosti v nehomogenih poroznih materialih, • določanje gostote dreves določene vrste v gozdu, • določanje lastnosti tal na izbranem področju, • analiza padavin na izbranem področju, • določanje meteoroloških parametrov v atmosferi, • določanje koncentracije onesnaževanja na onesnaženem področju, • arheološka izkopavanja, • seizmološke študije in • mnoge druge. 1.2 Podatki za prostorsko statistiko Osnovna razlika med podatki za običajno statistiko in prostorsko statistiko je v tem, da pri prostorski statistiki vsakemu elementu vzorca pripada informacija o legi v prostoru in/ali času, medtem ko pri običajni statistiki tega podatka ne upoštevamo. Na preprost način lahko to prikažemo v dveh vzporednih preglednicah: PREGLEDNICA 1.1: Primera podatkov pri prostorski in običajni statistiki Prostorska Običajna statistika statistika T [◦C] x[m] y[m] T [◦C] 22 112 70 22 15 142 107 15 27 47 78 27 18 30 145 18 · · · · · · · · · 1.3 Vrste obravnavnih slučajnih spremenljivk 3 V preglednicah 1.1 prikazujemo temperaturo v določenem času na obravnavanem območju. V obeh primerih obravnavamo isto slučajno spremenljivko, torej temperaturo. Razlika je le v tem, da so pri prostorski statistiki na voljo tudi informacije o tem, kje je bila temperatura merjena. 1.3 Vrste obravnavnih slučajnih spremenljivk Slučajne spremenljivke, ki jih obravnavamo v okviru prostorske statistike, so lahko diskretne ali zvezne, podobno kot smo to obravnavali v okviru osnov verjetnostnega računa in statistike. Informacija o legi v prostoru in času ima lahko različne značilnosti: • prostorska informacija je urejena ali neurejena (glej sliko 1.1), • prostorska informacija je lahko realizacija diskretne ali zvezne slučajne spremenljivke, • prostorska informacija je lahko podana v točkah (na primer vodnjak ali vrtina), na krivuljah (na primer cesta ali vodotok), ravninskih območjih (na primer parcela, občina ali država) ali po prostorskih območjih (na primer nahajališče rude). Prostorsko urejeni podatki Prostorsko neurejeni podatki - enakomerno Prostorsko neurejeni podatki - v gručah SLIKA 1.1: Različne razporeditve točk, v katerih so merjeni podatki Vrste prostorskih podatkov so: • geostatistični podatki (ang. geostatistical data), • mrežna struktura podatkov (ang. lattice data), • prostorski vzorci točk (ang. point patterns). 4 1 Uvod 1.3.1 Geostatistični podatki Podatki so pripravljeni tako, da so vrednosti slučajne spremenljivke, ki jo obravnavamo, merjene v točkah v prostoru, ki so lahko razporejene urejeno ali neurejeno. Tu opisujemo nekaj primerov geostatističnih podatkov, ki jih bralec lahko uporabi pri svojih analizah. Podatke lahko dobite na spletnih straneh UL FGG. • Podatki KIS - podatki o pedoloških profilih predstavljajo vrednosti treh količin na območju Slovenije na globini 0 - 20 cm: vrednost pH (v KCl), količina izmenljivega fosforja P (mg/100 g) in količina izmenljivega kalija K (mg/100 g) v tleh. Podatke so pripravili raziskovalci Kmetijskega inštituta Slovenije (B. Vrščaj in J. Bergant). Na sliki 1.2 so z rdečo označene točke, kjer količina prostega fosforja presega 25 mg/100 g. Lege merskih točk v Sloveniji 200 000 150 000 y [m] 100 000 50 000 400 000 450 000 500 000 550 000 600 000 x [m] SLIKA 1.2: Razpored točk na območju Slovenije z označenimi točkami z višjo koncentracijo fosforja • Podatki GASA - Geostatistical Association of Southern Africa so pridobljeni z meritvami vseb-nosti zlata v rudnini (g/t) in debeline plasti rudnine (cm) v globokih vrtinah v Južnoafriški repub-liki. Meritve so predstavljale osnovo za načrtovanje rudnika zlata Wits. • Podatki Wolfcamp predstavljajo meritve vodnega tlaka v 85 vodnjakih v New Mexicu in Texasu. Meritve je naročil Office for Nuclear Waste Isolation ob načrtovanju odlagališča radioaktivnih odpadkov. Meritve so uporabili za oceno, kako hitro bi morebitno puščanje radioaktivnih odpadkov onesnažilo vodo na tem območju. • Podatke Dioxin so izmerili na območju okoli ceste v Missouriju, kjer je iz tovornjaka iztekla 1.3 Vrste obravnavnih slučajnih spremenljivk 5 neznana količina strupene tekočine. V želji po racionalizaciji meritev so pripravili načrt meritev v urejenem vzorcu okoli ceste. Merili so koncentracijo strupene snovi TCDD (tetraklordibenzo-p-dioksin) v µg/kg. • Podatki Scallops predstavljajo podatke o številu ulovljenih školjk na priobalnem območju vzhodne obale ZDA. Poleg koordinat, kjer je bil vzet vzorec, so prešteli število školjk, manjših in večjih od 7 cm, ter skupno število školjk. • Podatke Jura so zbrali na Swiss Federal Institute of Technology v Lozani (fr. Lousanne) in predstavljajo značilne geostatistične podatke. Na znanih koordinatah 359 točk v švicarskem kantonu Jura so v tleh merili koncentracije kadmija, kobalta, kroma, bakra, niklja, svinca in cinka. V teh točkah so določili tudi tip kamnine (pet različnih kamnin) in trenutno rabo tal (gozd, pašnik, travnik ali njiva), glej sliko 1.3. Podatke o pokrajini Jura bomo večkrat uporabili pri obravnavanju različnih geostatističnih metod. Različne rabe zemljišč na območju kantona Jura Različni tipi kamnine na območju kanotna Jura gozd Argovian pašnik Kimmeridgian travnik Sequanian njiva Portlandian Quartarnary SLIKA 1.3: Razpored točk na območju kantona Jura z oznakami rabe zemljišča in tipa kamnine 6 1 Uvod 1.3.2 Podatki v mreži Mreže so lahko geometrijsko oziroma prostorsko pravilnih oblik, na primer kvadratne ali trikotne, ali pa nepravilnih oblik, kot so na primer geografska območja, občine, države in podobno. Prostorska informacija je podana topološko. To pomeni, da za posamezne elemente vemo, kateri drugi elementi so obravnavanemu elementu sosed. Na ta način lahko ugotavljamo, v kolikšni meri elementi vplivajo na svoje sosede. Podatki so pripravljeni tako, da imamo za vsak element seznam sosedov. Soseda opredelimo kot element, ki ima z obravnavanim skupno mejo. Med sosede lahko vključimo tudi tiste elemente, ki imajo skupno mejo z neposrednim sosedom. Na primer: Avstrija ima s Slovenijo skupno mejo in je po prvem kriteriju Sloveniji soseda. Če uporabimo še drugi kriterij, lahko med sosede štejemo tudi države, ki mejijo na Italijo, Avstrijo, Madžarsko ali Hrvaško, to so na primer Francija, Slovaška, Češka, Nemčija in Srbija (glej sliko 1.4). SLIKA 1.4: Slovenija s svojimi sosedami v prvem (modre) in drugem kolenu (zelene) 1.3 Vrste obravnavnih slučajnih spremenljivk 7 Sosede lahko definiramo tudi kot elemente, katerih koordinate, na primer geometrijsko središče elementa, so dovolj blizu obravnavanega elementa oziroma njegovega geometrijskega središča. Tak primer je prikazan na sliki 1.5, kjer smo med sosede Slovenije šteli Italijo, San Marino, Avstrijo, Češko, Madžarsko, Hrvaško ter Bosno in Hercegovino. SLIKA 1.5: Slovenija s svojimi sosedami glede na oddaljenost geometrijskih središč Tak način podajanja in obravnavanja prostorskih podatkov lahko uporabimo na primer pri obravnavanju pojavov, kot sta širjenje nalezljivih bolezni ali razvoj kriminalnih aktivnosti. 8 1 Uvod 1.3.3 Prostorski točkovni vzorci O prostorskih vzorcih govorimo, ko ima lega oziroma razporeditev točk v prostoru pomembno vlogo in jo želimo obravnavati. V tem primeru je eno prvih vprašanj, kako so točke razporejene, slučajno, v gručah ali urejeno, kar je grafično prikazano na slikah 1.1 in 1.6. Primere take obravnave najdemo v gozdarskih raziskavah, kjer na primer ugotavljajo razporejenost in lastnosti dreves smreke v mešanem gozdu, v bioloških raziskavah rasti celic v pogojih in-vitro in v mnogih drugih. SLIKA 1.6: Različni vzorci geometrijske porazdelitve dreves: zgoraj slučajno enakomerno, spodaj levo v gručah, spodaj desno urejen vzorec (foto: G. Turk) 2 Slučajne funkcije, polja in procesi V tem poglavju bomo obravnavali osnovne definicije slučajne funkcije, slučajnega polja oziroma procesa. Najprej pa bomo na kratko ponovili definicije, povezane s slučajno spremenljivko, tako da bodo poseb-nosti, ki so značilne za slučajno funkcijo, ki je osnova prostorske statistike, bolj razumljive. 2.1 Slučajne spremenljivke Slučajno spremenljivko označimo z veliko tiskano črko (na primer X, Y, Z), vrednosti, ki jih slučajna spremenljivka lahko zavzame, pa z malimi črkami (na primer x, y, z). Dogodke, povezane s slučajno spremenljivko, lahko zapišemo na primer takole: Z ≤ z, kar bi prebrali kot ,,slučajna spremenljivka Z zavzame vrednosti manjše ali enake od vrednosti z“, ali pa Z = 4, kar bi prebrali kot ,,slučajna spremenljivka zavzame vrednost 4“. Pomemben pojem pri obravnavanju slučajnih spremenljivk je zaloga vrednosti, ki predstavlja tisto množico vrednosti, ki jih slučajna spremenljivka lahko zavzame. Zaloga vrednosti je lahko diskretna množica vrednosti, v tem primeru govorimo o diskretni slučajni spremenljivki. Zaloga vrednosti je lahko tudi množica vrednosti na intervalu ali uniji intervalov, v tem primeru govorimo o zvezni slučajni spremenljivki. Primer zaloge vrednosti diskretne slučajne spremenljivke so naravna števila od 1 do 6 (z ∈ N, 1 ≤ z ≤ 6), kar lahko predstavlja zalogo vrednosti rezultata meta kocke. Primer zaloge vrednosti zvezne slučajne spremenljivke so pozitivna realna števila (z ∈ R, z > 0), kar lahko predstavlja čas do naslednjega potresa v Ljubljani. Za popolno predstavitev slučajne spremenljivke moramo poznati njeno porazdelitev. Le-ta je lahko podana na različne načine. S porazdelitveno funkcijo, ki je definirana z enačbo: FZ(z) = P [Z ≤ z], (2.1) lahko opišemo porazdelitev poljubne diskretne in zvezne slučajne spremenljivke. Če preberemo enačbo (2.1), je porazdelitvena funkcija FZ(z) slučajne spremenljivke Z definirana kot verjetnost, da slučajna spremenljivka Z zavzame vrednosti, ki so manjše ali enake vrednosti z. Pomembni lastnosti porazdelitvene funkcije, ki izvirata iz aksiomov verjetnostnega računa, sta, da je omejena med 0 in 1 ter da je nepadajoča: 10 2 Slučajne funkcije, polja in procesi 0 ≤ FZ(z) ≤ 1, lim FZ(z) = 0, z→−∞ lim FZ(z) = 1. z→∞ Porazdelitev diskretne slučajne spremenljivke lahko opišemo z verjetnostno funkcijo, ki je definirana z enačbo: pZ(zi) = P [Z = zi]. (2.2) Verjetnostna funkcija pZ(zi) je torej definirana kot verjetnost, da slučajna spremenljivka Z zavzame vrednost zi. Lastnosti verjetnostne funkcije povzamemo z naslednjima izrazoma: 0 ≤ pZ(zi) ≤ 1, n X pZ(zi) = 1, i=1 kjer je n število vrednosti v zalogi vrednosti slučajne spremenljivke Z. Porazdelitev zvezne slučajne spremenljivke lahko opišemo tudi z gostoto verjetnosti fZ(z), s katero lahko izračunamo verjetnost, da slučajna spremenljivka zavzame vrednost na določenem intervalu: Z b P [a < Z ≤ b] = fZ(z)dz. (2.3) a Lastnosti gostote verjetnosti lahko povzamemo z naslednjima izrazoma: 0 ≤ fZ(z), Z ∞ fZ(z)dz = 1. −∞ Pomembni sta tudi enačbi, ki povezujeta gostoto verjetnosti in porazdelitveno funkcijo Z z dFZ(z) FZ(z) = fZ(˜ z) d˜ z ←→ fZ(z) = . −∞ dz 2.2 Slučajni vektorji Če obravnavamo pojav, ki ga moramo opisati z več kot eno spremenljivko, le-te sestavimo v slučajni vektor. Slučajni vektor več slučajnih spremenljivk zapišemo takole: Z = {Z1, Z2, . . . , Zn}, 2.2 Slučajni vektorji 11 kjer je n število slučajnih spremenljivk v slučajnem vektorju oziroma število dimenzij, ki jih slučajni vektor opisuje. Podobno kot pri slučajni spremenljivki zaloga vrednosti slučajnega vektorja predstavlja vse vrednosti, ki jih slučajna spremenljivka lahko zavzame. Slučajne spremenljivke, ki sestavljajo slučajni vektor, so lahko tako zvezne kot diskretne. Tudi slučajni vektor v celoti predstavimo z njegovo porazdelitvijo. Porazdelitvena funkcija slučajnega vektorja je definirana z enačbo: FZ (z 1,Z2,...,Zn 1, z2, . . . , zn) = P [Z1 ≤ z1 ∩ Z2 ≤ z2 ∩ . . . ∩ Zn ≤ zn]. (2.4) Nekaj pomembnih lastnosti porazdelitvene funkcije slučajnega vektorja povzamemo z izrazi: 0 ≤ FZ (z 1,Z2,...,Zn 1, z2, . . . , zn) ≤ 1, lim FZ (z1, z2, . . . , zn) = 0 za vsak i = 1, . . . , n, z 1,Z2,...,Zn i→−∞ lim FZ (z1, z2, . . . , zn) = 1, z 1,Z2,...,Zn 1→∞ ∩ z2→∞ ∩... ∩ zn→∞ lim FZ (z1, z2, . . . , zn) = FZ (z2, z3, . . . , zn). z 1,Z2,...,Zn 2,Z3,...,Zn 1→∞ Z zadnjim izrazom definiramo porazdelitev dela slučajnega vektorja, taki porazdelitveni funkciji pravimo robna porazdelitvena funkcija. Podobno kot pri slučajni spremenljivki definiramo tudi verjetnostno funkcijo in gostoto verjetnosti slučajnega vektorja. Verjetnostna funkcija je: pZ (z 1,Z2,...,Zn 1i, z2j , . . . , znk) = P [Z1 = z1i ∩ Z2 = z2j ∩ . . . ∩ Zn = znk], (2.5) z gostoto verjetnosti pa lahko izračunamo verjetnost, da slučajni vektor leži na določenem intervalu: P [a1 < Z1 ≤ b1 ∩ a2 < Z2 ≤ b2 ∩ . . . ∩ an < Zn ≤ bn] = Z b1 Z b2 Z bn (2.6) = . . . fZ dz 1,Z2,...,Zn 1 dz2 . . . dzn. a1 a2 an Pri slučajnem vektorju lahko definiramo tudi pogojno porazdelitev, ki predstavlja porazdelitev slučajnega vektorja ob pogoju, da so za nekatere slučajne spremenljivke v slučajnem vektorju izpolnjeni določeni pogoji. Na primer, pogojno porazdelitveno funkcijo definiramo kot: FZ1|(Z2≤z2 ∩...∩ Zn≤zn)(z1) = P [(Z1 ≤ z1)|(Z2 ≤ z2 ∩ . . . ∩ Zn ≤ zn)] = P [Z1 ≤ z1 ∩ Z2 ≤ z2 ∩ . . . ∩ Zn ≤ zn] FZ (z1, z2, . . . , zn) (2.7) = = 1,Z2,...,Zn . P [Z2 ≤ z2 ∩ . . . ∩ Zn ≤ zn] FZ (z 2,...,Zn 2, . . . , zn) V primeru, da je slučajna spremenljivka Z1 neodvisna od drugih slučajnih spremenljivk v slučajnem vektorju, lahko zapišemo: FZ (z 1|(Z2≤z2 ∩...∩ Zn≤zn)(z1) = FZ1 1) −→ (2.8) −→ FZ (z (z (z 1,Z2,...,Zn 1, z2, . . . , zn) = FZ1 1)FZ2,...,Zn 2, . . . , zn), 12 2 Slučajne funkcije, polja in procesi kar pomeni, da lahko porazdelitveno funkcijo slučajnega vektorja, sestavljenega iz neodvisnih slučajnih spremenljivk, zapišemo kot produkt robnih porazdelitvenih funkcij. Torej, če so vse slučajne spremenljivke medsebojno neodvisne, lahko porazdelitveno funkcijo slučajnega vektorja zapišemo kot produkt vseh porazdelitvenih funkcij slučajnih spremenljivk, ki sestavljajo vektor: FZ (z (z (z (z 1,Z2,...,Zn 1, z2, . . . , zn) = FZ1 1)FZ2 2) . . . FZn n). (2.9) 2.3 Momenti slučajnih spremenljivk in vektorjev Momenti in centralni momenti slučajnih spremenljivk opisujejo lastnosti slučajnih spremenljivk in vektorjev. Najpomembnejši momenti slučajnih spremenljivk in vektorjev so: pričakovana vrednost, varianca in kovarianca ter simetričnost in sploščenost. V nadaljevanju bomo obravnavali le momente zveznih slučajnih spremenljivk in vektorjev. Pričakovana vrednost slučajne spremenljivke Z je moment prvega reda, definiran z izrazom: Z ∞ mZ = E[Z] = z fZ(z) dz, (2.10) −∞ in predstavlja osnovno informacijo o vrednosti slučajne spremenljivke. S črko E smo označili pričakovano vrednost, povzeto po angleških virih (angl. expected value). Varianca slučajne spremenljivke Z je centralni moment drugega reda, definiran z izrazom: Z ∞ σ2Z = var[Z] = (z − mZ)2 fZ(z) dz = E[Z2] − E[Z]2 (2.11) −∞ in predstavlja osnovno informacijo o razpršenosti slučajne spremenljivke. Centralna momenta tretje in četrte stopnje predstavljata meri za simetričnost in sploščenost: Z ∞ (3) µ = (z − m Z Z )3 fZ (z) dz, mera simetričnosti, −∞ (2.12) Z ∞ (4) µ = (z − m Z Z )4 fZ (z) dz, mera sploščenosti. −∞ Brezdimenzijska koeficienta simetričnosti in sploščenosti definiramo z enačbama: (3) µ γ Z 1Z = , koeficient simetričnosti, σ3Z (2.13) (4) µ γ Z 2Z = , kurtosis. σ4Z 2.4 Slučajne funkcije 13 Kovarianca med slučajnima spremenljivkama Z1 in Z2 je centralni moment drugega reda, definiran z izrazom: Z ∞ Z ∞ σZ = cov[Z (z )(z ) f (z 1,Z2 1, Z2] = 1 − mZ1 2 − mZ2 Z1,Z2 1, z2) dz1 dz2 −∞ −∞ (2.14) = E[Z1Z2] − E[Z1]E[Z2], in predstavlja mero za linearno povezanost med dvema slučajnima spremenljivkama. Brezdimenzijska mera linearne povezanosti je korelacijski koeficient ρZ , ki je definiran z enačbo: 1,Z2 σZ ρ 1,Z2 Z = (2.15) 1,Z2 σZ σ 1 Z2 in lahko zavzame vrednosti med −1 in 1: −1 ≤ ρZ ≤ 1. 1,Z2 2.4 Slučajne funkcije Slučajna funkcija je slučajna spremenljivka ali slučajni vektor, ki je odvisna tudi od lege v prostoru in časa. Če je slučajna funkcija odvisna le od lege v prostoru, ji pravimo slučajno polje. Če je slučajna funkcija odvisna tudi od časa, ji rečemo slučajni proces. Omejili se bomo le na zvezna slučajna polja in procese, saj le-ti predstavljajo pomembnejši del prostorske statistike. Skalarno slučajno polje zapišemo kot: Z(s), s ∈ n R , (2.16) skalarni proces pa kot: Z(s, t), s ∈ n R , t ∈ R, t > 0, oziroma Z(t), t ∈ R, t > 0, (2.17) kjer s predstavlja krajevni vektor do točke, kjer zapišemo obravnavano slučajno polje, t pa čas, pri katerem obravnavamo slučajni proces. Definirali smo torej slučajno spremenljivko Z(s) oziroma Z(s, t) ali Z(t), ki jo lahko obravnavamo podobno kot običajno slučajno spremenljivko. Pomen zaloge vrednosti je povsem enak kot pri slučajni spremenljivki, to so vrednosti, ki jih lahko ima slučajno polje. Primer slučajnega polja so trenutne temperature v Sloveniji, kjer lahko ugotovimo, da so tako pričakovane vrednosti kot tudi varianca in porazdelitev različni v različnih točkah. Najprej bomo obravnavali le skalarna slučajna polja, torej bomo obravnavali le slučajno spremenljivko, odvisno od lege. Porazdelitvena funkcija slučajnega polja FZ(s)(z) je: FZ(s)(z) = P [Z(s) ≤ z] (2.18) in pomeni verjetnost, da slučajno polje v točki, ki je določena s krajevnim vektorjem s, zavzame vrednost, manjšo ali enako z. Podobno kot moramo pri slučajnem vektorju {Z1, Z2} poznati porazdelitveno 14 2 Slučajne funkcije, polja in procesi funkcijo FZ (z (z (z 1,Z2 1, z2) in ne le robnih porazdelitvenih funkcij FZ1 1) in FZ2 2), moramo za celovit opis skalarnega slučajnega polja definirati porazdelitveno funkcijo slučajnega polja Z(s) za dve različni točki s1 in s2: FZ(s1),Z(s2)(z1, z2) = P [Z(s1) ≤ z1 ∩ Z(s2) ≤ z2], (2.19) kar je povsem analogno definiciji porazdelitvene funkcije slučajnega vektorja po enačbi (2.4). V primeru slučajnega vektorja smo porazdelitveno funkcijo definirali za dve ali več slučajnih spremenljivk, pri slučajnem polju pa smo v (2.19) definirali porazdelitveno funkcijo za isto slučajno spremenljivko v dveh različnih točkah. Gostoto verjetnosti slučajnega polja definiramo kot odvod porazdelitvene funkcije: dFZ(s)(z) fZ(s)(z) = (2.20) dz oziroma: ∂2FZ(s f 1),Z (s2)(z1, z2) Z(s . (2.21) 1),Z (s2)(z1, z2) = ∂z1∂z2 Podobno kot za slučajne spremenljivke lahko določimo še nekaj osnovnih momentov, s katerimi opišemo osnovne lastnosti slučajnega polja. Srednja vrednost slučajnega polja je: Z ∞ mZ(s) = E[Z(s)] = z fZ(s)(z) dz, (2.22) −∞ varianca slučajnega polja pa je: Z ∞ σ2 = var[Z(s)] = z − m 2 f Z(s) Z(s) Z(s)(z) dz = E[Z (s)2] − E[Z (s)]2 = C (s). (2.23) −∞ V primeru slučajnega polja lahko določimo še funkcije, ki jih pri slučajnih spremenljivkah nismo imeli. To so kovariančna funkcija, korelacijska funkcija in variogram. Kovariančna funkcija predstavlja kovarianco med slučajnim poljem v dveh različnih točkah s1 in s2: σZ(s1),Z(s2) = cov[Z(s1), Z(s2)] = Z ∞ Z ∞ = z 1 − m (2.24) Z(s z f 1) 2 − mZ(s2) Z(s1),Z(s2)(z1, z2) dz1 dz2 = −∞ −∞ = E[Z(s1)Z(s2)] − E[Z(s1)] E[Z(s2)] = C(s1, s2). V zadnjih dveh enačbah smo vpeljali oznaki za varianco C(s) in kovariančno funkcijo C(s1, s2), ki sta običajni pri prostorski statistiki. Korelacijsko funkcijo določimo na enak način kot korelacijski koeficient pri slučajnih spremenljivkah: cov[Z(s1), Z(s2)] C(s1, s2) ρ(s1, s2) = = . (2.25) pvar[Z(s p 1)] var[Z (s2)] C(s1) C(s2) 2.4 Slučajne funkcije 15 Variogram oziroma semivariogram je definiran kot polovica variance razlike vrednosti slučajnega polja v dveh različnih točkah. Zaradi take definicije (polovica variance) sta se uveljavili obe poimeno-vanji, a je prvo bolj pogosto in bo dosledno uporabljeno v nadaljevanju učbenika. h 2γ(s1, s2) = var[Z(s1)−Z(s2)] = E Z(s1)−Z(s2)2i−EZ(s1)−Z(s2)2. (2.26) Če se lastnosti slučajnega polja po prostoru ne spreminjajo, pomeni, da obravnavamo stacionarno slučajno polje. Primera stacionarnega in nestacionarnega slučajnega polja prikazujemo na sliki 2.1, kjer s koordinatama x, y opišemo lego točke, z pa opisuje vrednost slučajnega polja. z z x y x y SLIKA 2.1: Stacionarno (levo) in nestacionarno (desno) slučajno polje Stacionarnost slučajnega polja opredelimo s pogojem, da je porazdelitvena funkcija za eno ali dve točki neodvisna od lege točke, določene s krajevnim vektorjem s: FZ(s)(z) = P [Z(s) ≤ z] = FZ(z), (2.27) FZ(s1),Z(s2)(z1, z2) = P [Z(s) ≤ z1 ∩ Z(s + h) ≤ z2] = FZ(h)(z1, z2). Če veljata pogoja (2.27), sta od lege neodvisni tudi srednja vrednost in varianca, medtem ko sta kovariančna funkcija in variogram odvisni le od vektorja h, ki predstavlja vektor med dvema točkama (glej sliko 2.2): s2 = s1 + h, h = s2 − s1. (2.28) Zapišemo enačbe za srednjo vrednost, varianco, kovariančno funkcijo in variagram za stacionarna polja: mZ(s) = E[Z(s)] = mZ, σ2 = var[Z(s)] = σ2 Z(s) Z = C (0), C(s1, s2) = C(s1, s1 + h) = C(s, s + h) = C(h), (2.29) γ(s1, s2) = γ(s1, s1 + h) = γ(s, s + h) = γ(h), C(h) C(h) ρ(h) = = . pC(0)C(0) C(0) 16 2 Slučajne funkcije, polja in procesi 1.0 0.8 h y s1 s 0.6 2 0.4 0.2 s 0.00.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 x SLIKA 2.2: Krajevni vektorji do nekaterih točk v prostoru Varianco σ2z smo zapisali kot C(0), saj s primerjavo enačb (2.22) in (2.24) vidimo, da lahko varianco zapišemo kot kovariančno funkcijo, zapisano za isto točko: σ2 = E[Z(s)Z(s)] − E[Z(s)]E[Z(s)] = C(s, s) = C(s, s + 0) = C(0). (2.30) Z(s) Za stacionarno polje lahko izpeljemo tudi zvezo med variogramom in kovariančno funkcijo. Zapišemo enačbo variograma, izraženo s pričakovanimi vrednostmi: h 2 γ(h) = E Z(s) − Z(s + h)2i − EZ(s) − Z(s + h)2 = (v desnem izrazu upoštevamo linearnost operatorja E[·]) = EZ(s)2 − 2 Z(s)Z(s + h) + Z(s + h)2 − E[Z(s)] − E[Z(s + h)]2 = (v levem izrazu upoštevamo linearnost operatorja E[·], zaradi stacionarnosti polja je E[Z(s)] = E[Z(s + h)] = mZ, zato se zadnja dva člena odštejeta) = EZ(s)2 − m2 − Z 2 E[Z(s)Z(s + h)] + 2 m2Z + EZ(s + h)2 − m2Z = (prišteli in odšteli smo 2 m2Z, da bomo lahko tvorili varianco in kovariančno funkcijo) = C(0) − 2 C(h) + C(0) = 2 C(0) − 2 C(h). Sedaj lahko zapišemo končni rezultat: γ(h) = C(0) − C(h) in C(h) = C(0) − γ(h). (2.31) 2.4 Slučajne funkcije 17 Zapišemo lahko tudi zvezo med korelacijsko funkcijo in variogramom: C(h) C(0) − γ(h) γ(h) ρ(h) = = = 1 − . (2.32) C(0) C(0) C(0) Za običajna slučajna polja velja, da kovariančna funkcija pada z naraščanjem razdalje med točkami, obratno pa variogram z naraščanjem razdalje med točkama običajno narašča in se približuje vrednosti C(0). Zapišemo lahko: lim C(h) = 0 in lim γ(h) = C(0). (2.33) |h|→∞ |h|→∞ Nekaj značilnih kovariančnih funkcij in variogramov prikazujemo na sliki 2.3. Kovariančna funkcija Variogram močno prostorsko povezano srednje prostorsko povezano C(h) prostorsko nepovezano γ(h) močno prostorsko povezano srednje prostorsko povezano prostorsko nepovezano h h SLIKA 2.3: Kovariančna funkcija in variogram V grafu variograma na sliki 2.4 lahko prepoznamo nekaj značilnih lastnosti. Variogram običajno narašča, saj pričakujemo, da bo varianca razlike med vrednostima slučajnega polja z razdaljo naraščala do neke meje, ki jo imenujemo prag (ang. sill). Opazimo lahko tudi, da je prag dosežen na neki določeni razdalji, do katere lahko prepoznamo prostorsko odvisnost. Tej razdalji rečemo doseg (ang. range) slučajnega polja in predstavlja razdaljo med točkama, kjer vrednost v eni točki še opazno vpliva na vrednost druge točke. Zadnja in morda najbolj presenetljiva lastnost variograma je tako imenovan učinek zrna (ang. nugget effect), ki odraža variabilnost vrednosti slučajnega polja v točki. Pričakovali bi, da je varianca razlike med vrednostima slučajnega polja v isti točki enaka nič, a je zaradi merskih napak in možnih nezveznosti le-ta od nič različna. Ime tega učinka izvira iz rudarstva, ki je bilo eno od področij, v okviru katerega se je prostorska statistika hitro razvijala. Predstavljamo si zrno zlata, ki ima seveda povsem druge lastnosti kot kamnina, v kateri zrno leži. Če bi torej merili lastnosti na meji med zrnom in kamnino, bi morali ugotoviti nezveznost oziroma variabilnost vrednosti v isti točki. 18 2 Slučajne funkcije, polja in procesi Variogram doseg γ (h) prag učinek "zrna" h SLIKA 2.4: Variogram in njegove značilnosti 2.5 Vektorske slučajne funkcije Kot je omenjeno v razdelku Geostatistični podatki in bomo podrobneje spoznali v naslednjem poglavju, lahko podatki obsegajo vrednosti več slučajnih spremenljivk. Tako na primer geostatistični podatki s področja Jura v Švici, ki jih bomo uporabili v različnih analizah prostorske statistike, poleg opisnih podatkov, kot so raba tal in tip kamnine, vsebuje tudi numerične podatke o koncentraciji različnih kovin (Cd, Cu, Pb, Co, Cr, Ni in Zn) v izbranih točkah. To je zelo značilen primer slučajnega vektorskega polja oziroma slučajne vektorske funkcije. Vektorsko slučajno polje zapišemo kot: {Z N n 1(s), Z2(s), Z3(s), . . . , ZN (s)} ∈ R , s ∈ R , (2.34) vektorski slučajni proces pa kot: {Z N n 1(s, t), Z2(s, t), Z3(s, t), . . . , ZN (s, t)} ∈ R , s ∈ R , t ∈ R, t > 0, oziroma (2.35) {Z N 1(t), Z2(t), Z3(t), . . . , ZN (t)} ∈ R , t ∈ R, t > 0, kjer s predstavlja krajevni vektor v n-razsežnem prostoru do točke, kjer zapišemo obravnavano slučajno polje, t pa čas, pri katerem obravnavamo slučajni proces. Število N predstavlja razsežnost slučajnega vektorskega polja, torej število slučajnih spremenljivk, ki so vključene v slučajni vektor oziroma slučajno vektorsko funkcijo. V primeru podatkov iz Jure obravnavamo sedem različnih slučajnih spremenljivk (N = 7) v dvorazsežnem prostoru geografske dolžine in širine (n = 2). Zaradi preprostejšega zapisa bomo v nadaljevanju obravnavali le vektorska slučajna polja z dvema slučajnima spremenljivkama (N = 2), saj v takem polju lahko opišemo vse pomembne značilnosti 2.5 Vektorske slučajne funkcije 19 analize vektorskih polj in je posplošitev na večje število slučajnih spremenljivk relativno preprosta. Zapišimo vektorsko polje dveh slučajnih spremenljivk: {Z n 1(s), Z2(s)}, s ∈ R . (2.36) Ker obravnavamo zvezno vektorsko slučajno polje, je zaloga vrednosti neko zvezno območje na 2 R . Porazdelitvena funkcija vektorskega slučajnega polja FZ1(s),Z2(s)(z1, z2) je: FZ1(s),Z2(s)(z1, z2) = P [Z1(s) ≤ z1 ∩ Z2(s) ≤ z2] (2.37) in pomeni verjetnost, da slučajno vektorsko polje v točki s, zavzame vrednost na območju Z1 ≤ z1 ∩ Z2 ≤ z2. Podobno kot pri skalarnem slučajnem polju lahko v analizo vpeljemo dve točki, določeni z vektorjema s1 in s2, in definiramo novo porazdelitveno funkcijo: FZ1(s1),Z2(s2)(z1, z2) = P [Z1(s1) ≤ z1 ∩ Z2(s2) ≤ z2], (2.38) s čimer postavljamo osnovo za obravnavo dveh različnih slučajnih spremenljivk Z1 in Z2 v dveh različnih točkah s1 in s2. Gostoto verjetnosti vektorskega slučajnega polja definiramo kot odvod porazdelitvene funkcije: ∂2FZ f 1(s),Z2(s)(z1, z2) Z , (2.39) 1(s),Z2(s)(z1, z2) = ∂z1∂z2 oziroma ∂2FZ f 1(s1),Z2(s2)(z1, z2) Z . (2.40) 1(s1),Z2(s2)(z1, z2) = ∂z1∂z2 Podobno kot za skalarno slučajno polje lahko za vektorsko slučajno polje določimo osnovne momente, s katerimi opišemo osnovne lastnosti vektorskega slučajnega polja. Srednji vrednosti komponent vektorskega slučajnega polja sta Z ∞ Z ∞ mZ z 1(s) = E [Z1(s)] = 1 fZ1(s),Z2(s)(z1, z2) dz1 dz2, −∞ −∞ (2.41) Z ∞ Z ∞ mZ z 2(s) = E [Z2(s)] = 2 fZ1(s),Z2(s)(z1, z2) dz1 dz2, −∞ −∞ varianci komponent vektorskega slučajnega polja pa sta Z ∞ Z ∞ σ2 = var[Z z 2 f Z 1(s)] = 1 − m 1(s) Z1(s) Z1(s),Z2(s)(z1, z2) dz1 dz2, −∞ −∞ (2.42) Z ∞ Z ∞ σ2 = var[Z z 2 f Z 2(s)] = 2 − m 2(s) Z2(s) Z1(s),Z2(s)(z1, z2) dz1 dz2, −∞ −∞ 20 2 Slučajne funkcije, polja in procesi ki ju lahko zapišemo tudi takole: σ2 = E[Z Z 1(s)2] − E[Z1(s)]2 = C1(s), 1(s) (2.43) σ2 = E[Z Z 2(s)2] − E[Z2(s)]2 = C2(s). 2(s) Analogno kot pri skalarnem slučajnem polju določimo tudi kovariančni funkciji in variograma za kom-ponenti vektorskega slučajnega polja: σZ1(s1),Z1(s2) = cov[Z1(s1), Z1(s2)] = Z ∞ Z ∞ = z 1 − mZ z f (2.44) 1(s1) 2 − mZ1(s2) Z1(s1),Z1(s2)(z1, z2) dz1 dz2 = −∞ −∞ = E[Z1(s1)Z1(s2)] − E[Z1(s1)] E[Z1(s2)] = C1(s1, s2), σZ2(s1),Z2(s2) = cov[Z2(s1), Z2(s2)] = Z ∞ Z ∞ = z 1 − mZ z f (2.45) 2(s1) 2 − mZ2(s2) Z2(s1),Z2(s2)(z1, z2) dz1 dz2 = −∞ −∞ = E[Z2(s1)Z2(s2)] − E[Z2(s1)] E[Z2(s2)] = C2(s1, s2), h 2γ1(s1, s2) = var[Z1(s1) − Z1(s2)] = E Z1(s1) − Z1(s2)2i − EZ1(s1) − Z1(s2)2, (2.46) h 2γ2(s1, s2) = var[Z2(s1) − Z2(s2)] = E Z2(s1) − Z2(s2)2i − EZ2(s1) − Z2(s2)2. (2.47) Pri vektorskem slučajnem polju lahko definiramo še križno kovariančno funkcijo, ki je definirana kot kovarianca slučajne spremenljivke Z1 v točki s1 s slučajno spremenljivko Z2 v točki s2: σZ1(s1),Z2(s2) = cov[Z1(s1), Z2(s2)] = Z ∞ Z ∞ = z 1 − m (2.48) Z z f 1(s1) 2 − mZ2(s2) Z1(s1),Z2(s2)(z1, z2) dz1 dz2 = −∞ −∞ = E[Z1(s1)Z2(s2)] − E[Z1(s1)] E[Z2(s2)] = C12(s1, s2) in križni variogram oziroma kovariogram, ki je definiran kot kovarianca razlik (Z1(s1) − Z1(s2)) in (Z2(s1) − Z2(s2)): 2 γ12(s1, s2) = cov[(Z1(s1) − Z1(s2)), (Z2(s1) − Z2(s2))] = h i = E Z (2.49) 1(s1) − Z1(s2) Z2(s1) − Z2(s2) − − EZ1(s1) − Z1(s2) EZ2(s1) − Z2(s2). 2.5 Vektorske slučajne funkcije 21 Podobno kot pri skalarnih slučajnih poljih lahko tudi pri vektorskih slučajnih poljih opredelimo stacionarnost, za taka vektorska slučajna polja velja: mZ , m , 1(s) = E [Z1(s)] = mZ1 Z2(s) = E[Z2(s)] = mZ2 σ2 = var[Z = C = var[Z = C Z 1(s)] = σ2 1(0), σ2 2(s)] = σ2 2(0), 1(s) Z1 Z2(s) Z2 C1(s1, s2) = C1(s1, s1 + h) = C1(h), C2(s1, s2) = C2(s1, s1 + h) = C2(h) (2.50) γ1(s1, s2) = γ1(s1, s1 + h) = γ1(h), γ2(s1, s2) = γ2(s1, s1 + h) = γ2(h), C12(s1, s2) = C12(s1, s1 + h) = C12(h), γ12(s1, s2) = γ12(s1, s1 + h) = γ12(h). Z upoštevanjem enačbe (2.31) lahko zapišemo zveze med variogramoma γ1(h), γ2(h) in korelacijskima funkcijama C1(h), C2(h): γ1(h) = C1(0) − C1(h) in γ2(h) = C2(0) − C2(h). (2.51) Za stacionarno vektorsko polje lahko izpeljemo tudi zvezo med kovariogramom in križno kovariančno funkcijo. Zapišemo enačbo kovariograma, izraženo s pričakovanimi vrednostmi: h i 2 γ12(h) = E Z1(s) − Z1(s + h) Z2(s) − Z2(s + h) − − EZ1(s) − Z1(s + h)EZ2(s) − Z2(s + h) = (upoštevamo linearnost operatorja E[·]) = EZ1(s)Z2(s) − Z1(s)Z2(s + h) − Z1(s + h)Z2(s) + Z1(s + h)Z2(s + h)− − E[Z1(s)] − E[Z1(s + h)] E[Z2(s)] − E[Z2(s + h)] = (upoštevamo linearnost operatorja E[·], zaradi stacionarnosti polja je E[Z1(s)] = E[Z1(s + h)] in E[Z2(s)] = E[Z2(s + h)], zato je drugi del izraza enak nič) = EZ1(s)Z2(s) − mZ m − E[Z m − 1 Z2 1(s)Z2(s + h)] + mZ1 Z2 − E[Z1(s + h)Z2(s)] + mZ m + EZ m = 1 Z2 1(s + h)EZ2(s + h) − mZ1 Z2 (prišteli in odšteli smo 2 mZ m , 1 Z2 da bomo lahko tvorili kovariance in križno kovariančno funkcijo) = C12(0) − C12(h) − C21(h) + C12(0). Sedaj lahko zapišemo končni rezultat: 1 γ12(h) = C12(0) − C12(h) + C21(h). (2.52) 2 Običajno predpostavimo simetričnost C12(h) = C21(h), zato sledi enaka zveza kot pri zvezi med vario-grami in korelacijskimi funkcijami (2.51) γ12(h) = C12(0) − C12(h). (2.53) 22 2 Slučajne funkcije, polja in procesi 3 Geostatistika V prejšnjem poglavju smo opisali slučajno polje oziroma slučajni proces in njihove osnovne značilnosti, porazdelitve, momente, mere prostorske povezanosti. V tem poglavju prehajamo na statistično analizo, obravnavali bomo vzorce slučajnih polj oziroma procesov in opisali ustrezne statistične metode. 3.1 Geostatistični podatki Geostatistične podatke sestavljajo podatki o legi točke, kjer je bila opravljena meritev količine, ki jo obravnavamo, in vrednosti meritev te količine. Primer, ki ga bomo uporabili za prikaz različnih geostatističnih metod, so podatki iz severozahodnega švicarskega kantona Jura na meji s Francijo. Tam je Švicarski zvezni inštitut za tehnologijo iz Lozane (EPFL) na približno 14 km2 v 359 točkah ocenil tip kamnine, rabo tal in izmeril koncentracijo težkih kovin v tleh [mg/kg]. Geostatistični podatki so v tem primeru sestavljeni iz koordinat merskih točk, podatkov o tipu kamnine in rabi tal ter rezultatov meritev koncentracij težkih kovin v tleh. Lege merskih točk 6 Koncentracije Cu in Zn v prsti 140 ] 5 Baker km 120 [ 4 y 100 Cink 3 80 2 60 Frekvenca 40 1 20 0 0 0 1 2 3 4 5 0 50 100 150 x [km] Koncentracija [mg/kg] SLIKA 3.1: Razpored merskih točk in histogram koncentracije dveh kovin Na sliki 3.1 je prikazan razpored merskih točk. Vidimo, da je večina točk razporejenih v pravilni mreži, poleg tega pa je dodanih še nekaj merskih mest, ki so dodana v neurejeni obliki. Lahko predvidevamo, da so bile te dodatne merske točke postavljene v bližino potencialnih onesnaževalcev ali pa v bližino zajetja 24 3 Geostatistika pitne vode, bližino vrtcev ali podobno. Zanimiva je tudi primerjava v porazdelitvi koncentracije bakra in cinka, saj vidimo, da so koncentracije cinka bistveno višje, njihova porazdelitev pa je bolj simetrična od porazdelitve koncentracije bakra. Lastnosti slučajnih polj koncentracije bakra in cinka bomo prikazali v nadaljevanju. Začetnih nekaj vrstic podatkov prikazujemo v preglednici 3.1. Razlikovali so med petimi tipi kamnine: prvi štirje (Argovian, Kimmeridgian, Sequanian in Portlandian) pripadajo jurskim apnencem, peta (Quartarnary) pa je mlajša kamnina iz obdobja kvartarja. Rabo tal so klasificirali v štiri razrede: gozd, pašnik, travnik in polja. PREGLEDNICA 3.1: Prvih nekaj vrstic podatkov za območje Jure x y Tip Raba Cd Cu Pb Co Cr Ni Zn [km] [km] kamnine tal [mg/kg] [mg/kg] [mg/kg] [mg/kg] [mg/kg] [mg/kg] [mg/kg] 2.386 3.077 3 3 1.740 25.72 77.36 9.32 38.32 21.32 92.56 2.544 1.972 2 2 1.335 24.76 77.88 10.00 40.20 29.72 73.56 2.807 3.347 3 2 1.610 8.88 30.80 10.60 47.00 21.40 64.80 4.308 1.933 2 3 2.150 22.70 56.40 11.92 43.52 29.72 90.00 4.383 1.081 5 3 1.565 34.32 66.40 16.32 38.52 26.20 88.40 3.244 4.519 5 3 1.145 31.28 72.40 3.50 40.40 22.04 75.20 3.925 3.785 5 3 0.894 27.44 60.00 15.08 30.52 21.76 72.40 2.116 3.498 1 3 0.525 66.12 141.00 4.20 25.40 9.72 72.08 1.842 0.989 1 3 0.240 22.32 52.40 4.52 27.96 11.32 56.40 1.709 1.843 3 3 0.625 18.72 41.60 12.08 33.32 16.88 75.60 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · Iz podatkov lahko določimo ocene srednjih vrednosti (povprečja) in variance (vzorčne variance) ter nekaj drugih statističnih podatkov, kar prikazujemo v preglednici 3.2. Iz te lahko preberemo, da naj-višje povprečne koncentracije izkazujeta svinec in cink, ki imata tudi največjo razpršenost. Verjetnostne porazdelitve koncentracije kobalta, kroma in niklja so najbližje normalni porazdelitvi, saj so te porazdelitve precej simetrične (koeficient simetričnosti je blizu nič) po obliki pa niso niti špičaste niti sploščene (kurtosis je blizu vrednosti tri), kar je značilnost normalne porazdelitve. 3.2 Priprava podatkov 25 PREGLEDNICA 3.2: Nekaj osnovnih statističnih podatkov Slučajne spremenljivke - elementi Statistika Cd Cu Pb Co Cr Ni Zn Povprečje 1.30 23.6 54.6 9.48 35.0 20.0 75.9 Standardni odklon 0.873 22.3 33.1 3.63 10.7 8.14 30.8 Mediana 1.11 17.2 46.8 9.84 34.8 20.6 73.6 Simetričnost 1.49 3.01 3.34 −0.0705 0.295 0.104 1.49 Kurtosis 5.66 14.7 18.6 2.72 3.32 3.17 7.70 5. centil 0.32 5.72 26.3 3.55 18.5 6.08 35.1 95. centil 3.15 65.2 116. 14.6 52.8 32.6 136. Minimum 0.135 3.55 18.7 1.55 3.32 1.98 25.0 Maksimum 5.13 166. 300. 22.7 70.0 53.2 260. Pomembna je tudi informacija, v kolikšni meri so slučajne spremenljivke povezane. Mera za linearno povezanost je kovarianca oziroma korelacijski koeficient. V preglednici 3.3 so prikazani korelacijski koeficienti, ki prikazujejo linearno povezanost med spremenljivkami. Korelacijski koeficient, ki se po absolutni vrednosti približuje vrednosti 1, predstavlja močno linearno povezani slučajni spremenljivki. Vidimo lahko, da sta močno linearno povezani koncentraciji bakra in svinca, zelo šibko pa sta povezni koncentraciji svinca in kadmija. PREGLEDNICA 3.3: Korelacijski koeficienti med obravnavanimi slučajnimi spremenljivkami Cd Cu Pb Co Cr Ni Zn Cd 1. 0.139 0.221 0.259 0.566 0.475 0.601 Cu 0.139 1. 0.824 0.183 0.215 0.222 0.657 Pb 0.221 0.824 1. 0.148 0.256 0.271 0.665 Co 0.259 0.183 0.148 1. 0.472 0.703 0.441 Cr 0.566 0.215 0.256 0.472 1. 0.705 0.61 Ni 0.475 0.222 0.271 0.703 0.705 1. 0.591 Zn 0.601 0.657 0.665 0.441 0.61 0.591 1. 3.2 Priprava podatkov Ob predpostavki, da je polje stacionarno, lahko iz podatkov določimo ocene srednjih vrednosti in varianc za vse merjene količine. Pri določitvi ocene kovariančne funkcije in variograma pa bomo morali poiskati pare točk z enako ali približno enako medsebojno oddaljenostjo. Zato izberemo vektor h z dolžino h = |h| in smernim kotom α. Za ta vektor poiščemo vse tiste pare točk, za katere velja, da je medsebojna 26 3 Geostatistika razdalja znotraj intervala h ± ∆h in da je smerni kot med prvo in drugo točko znotraj intervala α ± ∆α. Število teh parov označimo z N(h). Na sliki 3.2 prikazujemo množico 359 točk in lahko ugotovimo, da za izbrani vektor oddaljenosti med točkama h točko 17 lahko povežemo s točko 45 (obe sta obarvani v modro), medtem ko točka 18 za izbrani vekor h nima para (obarvana je rdeče). Na desni strani slike 3.2 je prikazanih vseh N (h) = 29 parov točk s približno enako medsebojno oddaljenostjo in približno enakim smernim kotom med točkama. Vidimo, da je lahko ena točka povezana z več točkami. V prikazanem primeru smo iskali le tiste pare točk z medsebojno oddaljenostjo h ± ∆h in s smernim kotom daljice od prve do druge točke v intervalu α ± ∆α. Pogosto predpostavimo, da je polje izotropno, torej da se lastnosti v različnih smereh ne razlikujejo. Tedaj pare točk s podobno oddaljenostjo iščemo ne glede na smerni kot daljice med dvema točkama. Tak izbor parov za istih 359 točk je prikazan na sliki 3.3. Čeprav smo velikost intervala ∆h zmanjšali za 8-krat, je število parov sedaj bistveno večje, N (h) = 117. Tudi na sliki 3.3 so z modro označene točke, ki so v parih, z rdečo pa tiste, ki nimajo para. Najpogosteje bomo obravnavali primere, ki niso odvisni od smeri vektorja h, zato bomo v nadaljevanju v teh primerih namesto vektorja h pri zapisu N (h) upoštevali le dolžino vektorja h in zapisali N (|h|) = N (h.) 6 6 Jura − iskanje parov Jura − iskanje parov h = 0.8 km N( h) = 29 parov ∆ h = 0.09 km h = 0.8 km 5 ∆α α = 0.4 h 5 α ∆ h ∆ h = 0.04 km ∆ α = 0.1 18 α = 0.4 ∆α = 0.05 4 45 4 17 y y 3 3 2 2 1 1 0 0 0 1 2 3 4 50 1 2 3 4 5 SLIKA 3.2: Iskanje parov točk, ki so oddaljene za približno h x x 3.2 Priprava podatkov 27 6 Jura − iskanje parov N( h) = 117 parov h = 0.8 km 5 ∆ h = 0.005 km 4 y 3 2 1 00 1 2 3 4 5 SLIKA 3.3: Pari točk s približno x oddaljenostjo h ne glede na smer daljice V vsaki točki sk poznamo tudi vrednost Zk slučajnega polja Z(s). Zato lahko ob iskanju parov za nek izbran h za vsak par točk na razdalji h zapišemo par vrednosti (Z1,i, Z2,i), kjer je Z1,i vrednost slučajnega polja v prvi točki i-tega para, Z2,i pa vrednost polja v drugi točki tega para. Na primer, če smo na sliki 3.2 ugotovili, da točki 17 in 45 predstavljata par dveh točk na ustrezni oddaljenosti glede na izbrani vektor h, potem je vrednost Z1,1 = Z17 (vrednost Z v prvi točki prvega para je enaka vrednosti Z17), vrednost v drugi točki pa je Z2,1 = Z45 (vrednost Z v drugi točki prvega para je enaka vrednosti Z45). S tem smo podatke na novo uredili tako, da smo iz geostatističih podatkov {xk, yk, Zk}, k = 1, . . . , n izluščili pare vrednosti {Z1,i, Z2,i}, i = 1, . . . , N (h) za določen vektor h, kot je prikazano v preglednici 3.4. V primeru s slike 3.2 je bilo vseh točk n = 359, ustreznih parov točk za izbrani vektor z dolžino h = 0.8 km in smerjo α = 0.4 pa je N (h) = 29. Zaradi preglednosti bomo v nadaljevanju vrednost slučajnega polja za določen par i pisali brez vejic med indeksi Z1,i → Z1i in Z2,i → Z2i. 28 3 Geostatistika PREGLEDNICA 3.4: Preurejevanje podatkov za geostatistično analizo Geostatistični podatki Pari vrednosti x1 y1 Z1 Z11 Z21 x2 y2 Z2 Z12 Z22 · · · · · · · · · −→ · · · · · · xk yk Zk Z1i Z2i · · · · · · · · · · · · · · · xn yn Zn Z1N(h) Z2N(h) S tako urejenimi podatki lahko grafično ali numerično prikažemo prostorsko povezanost. Najpreprostejši grafični prikaz je razsevni graf, ki ga lahko narišemo za različne razdalje h. Za vsako razdaljo h lahko izračunamo tudi oceno kovariančne in korelacijske funkcije ( ˆ C(h) in ˆ ρ(h)) ter oceno variograma ˆ γ(h). 3.3 Razsevni graf Z razsevnim grafom narišemo pare točk {Z1i, Z2i} za izbrane razdalje h med točkami. Na sliki 3.4 prikazujemo razsevni graf za krom za dve različni razdalji h med točkami. Če so točke razsevnega grafa bližje simetrijski osi prvega kvadranta, pomeni, da je slučajno polje za izbrani h prostorsko bolj povezano (glej sliko 3.4 za h = 0.2 km). Za prostorsko neodvisno slučajno polje je značilno, da so točke razsevnega grafa razporejene po celotnem območju in niso zbrane v bližini simetrijske osi (glej sliko 3.4 za h = 1.6 km). 80 80 Razsevni graf za Cr Razsevni graf za Cr h = 0.2 km h = 1.6 km 60 60 Z2 Z2 40 40 20 20 0 0 0 20 40 60 80 0 20 40 60 80 Z1 Z1 SLIKA 3.4: Razsevni graf za Cr pri h = 0.2 km in h = 1.6 km 3.4 Kovariančna funkcija 29 3.4 Kovariančna funkcija Za določitev ocene kovariančne in korelacijske funkcije moramo najprej oceniti srednje vrednosti in variance slučajnega polja. Te določimo z naslednjimi enačbami: N (h) 1 X ˆ m1 = Z1i, N (h) i=1 N (h) 1 X ˆ m2 = Z2i, N (h) i=1 (3.1) N (h) N (h) 1 X 1 X ˆ σ2 − 1 = (Z1i − ˆ m1)2 = Z2 ˆ m2 N (h) N (h) 1i 1, i=1 i=1 N (h) N (h) 1 X 1 X ˆ σ2 − 2 = (Z2i − ˆ m2)2 = Z2 ˆ m2 N (h) N (h) 2i 2. i=1 i=1 Predpostavka o stacionarnih poljih sicer zahteva, da sta enaki srednji vrednosti in varianci m1 = m2 in σ2 = σ2 1 2 , a zaradi naključne narave slučajnih procesov ne moremo pričakovati, da bosta enaki tudi oceni srednjih vrednosti in varianc v prvih oziroma drugih točkah parov. Velja torej: ˆ m1 ̸= ˆ m2, ˆ σ2 ̸ 1 = ˆ σ22. V primeru, da sta oceni srednje vrednosti in variance za prve točke parov zelo različni, predpostavka o stacionarnosti ne velja in moramo slučajno polje obravnavati kot nestacionarno. Oceno kovariančne funkcije oziroma vzorčno kovariančno funkcijo izračunamo po naslednji enačbi: N (h) N (h) ˆ 1 X 1 X C(h) = (Z1i − ˆ m1)(Z2i − ˆ m2) = Z1iZ2i − ˆ m1 ˆ m2, (3.2) N (h) N (h) i=1 i=1 oceno korelacijske funkcije pa po enačbi: ˆ C(h) ˆ ρ(h) = , kjer je − 1 ≤ ˆ ρ(h) ≤ 1. (3.3) ˆ σ1 ˆ σ2 Ocene kovariančne in korelacijske funkcije za h = 0.1, . . . , 2.0 km prikazujemo na sliki 3.5. Vidimo lahko, da je povezava nekoliko bolj izrazita za razdalje do 0.4 km, nato pa vrednost funkcije pade in prostorska povezanost slučajnega polja ni več značilna. Če primerjamo sliko 3.5 s sliko 3.6, kjer sta prikazani kovariančna in korelacijska funkcija za prostorsko neodvisno polje z enakima srednjo vrednostjo in varianco kot prostorsko povezano polje, vidimo, da imata v primeru neodvisnih polj tako kovariančna kot korelacijska funkcija za vse vrednosti h nizke vrednosti, ki so od nič različne le zaradi naključne narave slučajnega polja. 30 3 Geostatistika 40 Kovariančna funkcija za Cr 1.0 Korelacijska funkcija za Cr 30 0.5 C(h)20 ρ(h) 0.0 10 0 - 0.5 - 10 - 1.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 h h SLIKA 3.5: Oceni kovariančne in korelacijske funkcije za Cr 40 Kovariančna funkcija za Cr 1.0 Korelacijska funkcija za Cr 30 Prostorsko neodvisno polje Prostorsko neodvisno polje 0.5 C(h)20 ρ(h) 0.0 10 0 - 0.5 - 10 - 1.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 h h SLIKA 3.6: Oceni kovariančne in korelacijske funkcije za Cr (prostorsko neodvisno polje) 3.5 Variogram Oceno variograma oziroma vzorčni variogram izračunamo po naslednji enačbi: N (h) 1 X 2ˆ γ(h) = (Z1i − Z2i − ( ˆ m1 − ˆ m2))2, N (h) i=1 (zaradi stacionarnosti polja vzamemo, da je ( ˆ m1 − ˆ m2) = 0) (3.4) N (h) 1 X 2ˆ γ(h) = (Z1i − Z2i)2. N (h) i=1 Variogram predstavlja polovico variance razlike med vrednostjo slučajnega polja v dveh točkah na oddaljenosti h. Pričakujemo lahko, da bodo razlike med vrednostmi polja pri večjih razdaljah večje, torej da bo variogram naraščajoča funkcija h. Če pa je polje prostorsko neodvisno, lahko pričakujemo, da bo variogram konstanten. Na sliki 3.7 vidimo, da variogram pri dejanskih podatkih o koncentraciji kroma na področju Jure narašča, medtem ko je pri podatkih, ki so bili generirani kot prostorsko neodvisno polje, variogram konstanten v okviru nezanesljivosti zaradi naključne narave slučajnih polj. 3.6 Določitev empirične funkcije variograma 31 120 120 γ (h) γ (h) 80 80 40 40 Variogram za Cr Variogram za Cr Prostorsko neodvisno polje 0 0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 h h SLIKA 3.7: Oceni variograma za Cr (prostorsko odvisno in neodvisno polje) 3.6 Določitev empirične funkcije variograma V naslednjem poglavju bomo obravnavali krigiranje, pri katerem bomo potrebovali vrednosti variograma oziroma kovariančne funkcije v odvisnosti od razdalje med točkami h. Zato bomo iz vzorčnih vrednosti variograma, ki so podane po točkah, določili približno funkcijsko zvezo. Naštejmo nekaj značilnih funkcij, ki so glede na običajno obliko variograma primerne za njegovo aproksimacijo: g1(h) = 1 konstanta ( 3h − h3 h ≤ a g 2a 2a3 2(h) = sferična 1 h > a (3.5) g3(h) = 1 − e−3h/a eksponentna g4(h) = 1 − e−3h2/a2 Gaussova 1.2 1.0 g(h)0.8 g(h) = 1 0.6 3 h 1 - h3 h ≤ a g(h) = 2 a 2 a3 0.4 1 h > a 0.2 g(h) = 1 - e- 3 h/a g(h) = 1 - e- 3 h2/ a2 0.00.0 0.5 1.0 1.5 2.0 h SLIKA 3.8: Nekaj funkcij, primernih za aproksimacijo variograma Empirično funkcijo variograma izračunamo s seštevanjem funkcij iz (3.5) po naslednji enačbi: g(h) = b1g1(h) + b2g2(h, a2) + b3g3(h, a3) + · · · , (3.6) 32 3 Geostatistika pri čemer neznane koeficiente b1, b2, a2, . . . določimo po metodi najmanjših kvadratov. Če se na primer odločimo za kombinacijo konstante in eksponentne funkcije, je empirična funkcija variograma: g(h) = b1 + b2(1 − e−3h/a2), (3.7) če pa izberemo kombinacijo konstante in Gaussove funkcije, je empirična funkcija: g(h) = b1 + b2(1 − e−3h2/a22), (3.8) pri čemer moramo določiti še parametre b1, b2 in a2 tako, da se empirična funkcija čim bolj prilega točkovnim vrednostim vzorčnega variograma. Pri metodi najmanjših kvadratov moramo poiskati take parametre empiričnih funkcij variograma, da je vsota kvadratov razlik med vrednostmi vzorčnega variograma in izbrano funkcijo čim manjša: nv nv X X 2 SS = (γi − g(hi))2 = γi − b1 − b2(1 − e−3hi/a2) za eksponentno funkcijo, (3.9) i=1 i=1 kjer nv predstavlja število točk vzorčnega variograma. Vrednosti hi izberemo glede na lege točk, saj tako za premajhne kot za prevelike vrednosti hi ne moremo najti ustreznega para, torej je N (h) = 0 in kovariančne funkcije oziroma variograma ne moremo izračunati. Za določitev empirične kovariančne funkcije uporabimo zvezo med kovariančno funkcijo in variogramom (2.31) in upoštevamo, da vsota vseh bi predstavlja največjo vrednost variograma, ki je enaka varianci slučajnega polja C(0). Empirično kovariančno funkcijo zapišemo z enačbo: C(h) = (b1 + b2 + · · · ) − (b1 g1(h) + b2 g2(h, a2) + b3 g3(h, a3) + · · · ). (3.10) Učinek zrna je v empirični funkciji variograma opisan s parametrom b1, ki predstavlja konstantni del variograma. Sedaj določimo empirični variogram za krom za razdalje med točkami od 0 do 2.5 km ob uporabi dveh različnih nastavkov po enačbah (3.7) in (3.8). Rezultat aproksimacije so vrednosti parametrov za obe enačbi, ki jih podajamo v preglednici 3.5, aproksimacijo variograma z vzorčnim variogramom pa prikazujemo na sliki 3.9. PREGLEDNICA 3.5: Parametri aproksimacije variograma za krom po enačbah (3.7) in (3.8) Funkcija nv b1 b2 a2 Vsota kvadratov SS (3.7) 49 36.5 77.3 0.604 5529 (3.8) 49 63.1 50.1 0.498 5496 Iz preglednice 3.5 in s slike 3.9 lahko razberemo, da sta prag in doseg obeh funkcij variograma podobni, saj je ta vrednost pri eksponentni funkciji 113.8, pri Gaussovi pa 113.1. Razlika v učinku zrna je bolj izrazita: v primeru eksponentne funkcije je enak 36.5, v primeru Gaussove pa 63.1. Vsota kvadratov odstopanj SS je v obeh primerih zelo podobna (5529 in 5496), Gaussova aproksimacija kaže nekoliko boljše prileganje, a je razlika zanemarljiva. 3.7 Geostatistična analiza z indikatorji 33 120 g(h) 80 40 g(h) - konstanta in eksponentna (3.8) g(h) - konstanta in Gaussova (3.9) 0 0.0 0.5 1.0 1.5 2.0 2.5 h SLIKA 3.9: Aproksimacija variograma 3.7 Geostatistična analiza z indikatorji Pri analizi geostatističnih podatkov nas pogosto zanimajo le vrednosti, ki so manjše ali večje od izbrane kritične vrednosti. Tako nas na primer lahko zanima, ali so vrednosti, ki presegajo zakonsko dovoljene koncentracije onesnaženja, prostorsko odvisne ali ne. Podobno nas lahko zanima, ali je vsebnost rude v kamnini, manjša od določene kritične vrednosti, prostorsko odvisna in kako. Če želimo opraviti tako analizo, moramo iz osnovnih geostatističnih podatkov določiti takoimenovane indikatorje po naslednji enačbi: 1 za Z i k ≤ zkrit, k = (3.11) 0 za Zk > zkrit, ali pa: 1 za Z i k ≥ zkrit, k = (3.12) 0 za Zk < zkrit. Geostatistične podatke smo sedaj razširili za indikatorje ik(zkrit) glede na izbrano kritično vrednost zkrit. Ob določitvi parov za izbrano razdaljo med točkami h dobimo nove pare vrednosti {i1i, i2i}, kot prikazujemo v preglednici 3.6. PREGLEDNICA 3.6: Preurejevanje podatkov za geostatistično analizo z indikatorji Geostatistični podatki Pari indikatorjev x1 y1 Z1 i1 i11 i21 x2 y2 Z2 i2 i12 i22 · · · · · · · · · · · · −→ · · · · · · xk yk Zk ik i1i i2i · · · · · · · · · · · · · · · · · · xn yn Zn in i1N(h) i2N(h) Glede na to, da indikatorji lahko dosežejo le dve vrednosti (ena ali nič), razsevnega grafa ne rišemo, saj 34 3 Geostatistika lahko vsebuje največ 4 točke ({0, 0}, {1, 0}, {0, 1}, {1, 1}). Zato pa lahko izračunamo vse druge mere: srednjo vrednost, varianco, kovariančno in korelacijsko funkcijo ter variogram: N (h) 1 X ˆ mI1 = i1i, N (h) i=1 (3.13) N (h) 1 X ˆ mI2 = i2i, N (h) i=1 N (h) 1 X ˆ σ2 − I1 = i2 ˆ m2 N (h) 1i I1, i=1 N (h) 1 X ˆ σ2 − I2 = i2 ˆ m2 N (h) 2i I2, i=1 N (h) ˆ 1 X (3.14) CI (h) = i1ii2i − ˆ mI1 ˆ mI2, N (h) i=1 ˆ CI (h) ˆ ρI (h) = , ˆ σI1 ˆ σI2 N (h) 1 X 2ˆ γI (h) = (i1i − i2i)2. N (h) i=1 Na slikah 3.10 in 3.11 prikazujemo ocene za kovariančno in korelacijsko funkcijo ter variograma za kritično vrednost zkrit = 25 mg/kg koncentracije kroma, pri čemer smo indikatorje določili po enačbi (3.12). Vidimo lahko, da je za izbrano kritično vrednost 25 mg/kg prostorska povezanost podobna kot, če obravnavamo vrednosti koncentracije. 0.10 1.0 Kovariančna funkcija Korelacijska funkcija 0.08 C 0.5 ρ I(h)0.06 CI (h) za Cr I (h) za Cr ρ I (h) 0.04 0.0 0.02 zkrit = 25 mg/kg 0.00 - 0.5 zkrit = 25 mg/kg - 0.020.0 0.5 1.0 1.5 2.0 - 1.0 h 0.0 0.5 1.0 1.5 2.0 h SLIKA 3.10: Oceni kovariančne in korelacijske funkcije z indikatorji za Cr 3.8 Preizkušanje domnev o prostorski povezanosti 35 0.20 0.15 γI(h) 0.10 z 0.05 krit = 25 mg/kg Variogram ℽI(h) za Cr 0.000.0 0.5 1.0 1.5 2.0 h SLIKA 3.11: Ocena variograma z indikatorji za Cr Na sliki 3.12 prikazujemo kovariančno funkcijo z indikatorji za različne kritične vrednost. Opazimo lahko, da je prostorska povezanost izrazitejša za kritične vrednosti med 25 in 35mg/kg, za manjše ali večje kritične vrednosti pa je povezanost manjša. 0.10 Kovariančna funkcija 0.08 CI(h)0.06 CI (h) za Cr 0.04 0.02 zkrit = 25 mg/kg 0.00 - 0.020.0 0.5 1.0 1.5 2.0 h SLIKA 3.12: Kovariančna funkcija z indikatorji za Cr za različne vrednosti zkrit 3.8 Preizkušanje domnev o prostorski povezanosti Geostatistične analize, ki smo jih opisali do sedaj, so pokazale prostorsko odvisnost preko oblike razsevnega grafa, preko oblike kovariančne, korelacijske funkcije ali variograma. Mnenje o tem, ali je slučajno polje prostorsko povezano ali ne, pa je prepuščeno subjektivni presoji grafičnih prikazov. V nadaljevanju bomo pokazali, da lahko pri presoji o prostorski povezanosti slučajnega polja uporabimo običajno statistično orodje: preizkušanje domnev (tudi testiranje hipotez, angl. hypothesis testing). Postopek preizkušanja domnev o prostorski povezanosti slučajnega polja je povsem enak kot pri vseh drugih preizkušanjih statističnih domnev: 36 3 Geostatistika 1) Postavitev ničelne H0 in alternativne domneve H1. Alternativna domneva ne sme biti združljiva z ničelno. 2) Izbira največjega tveganja α, ki ga dovolimo pri zavrnitvi ničelne domneve (običajne vrednosti so α = 0.01, α = 0.05 ali α = 0.10). 3) Zapis statistike, ki je skladna z ničelno domnevo. 4) Določitev mej kritičnega območja oziroma območja zavrnitve ničelne domneve. 5) Glede na izračunano statistiko v primerjavi z mejami kritičnega območja odločitev o zavrnitvi ničelne domneve in zapis končne izjave. Ničelna domneva je pri vseh preizkusih o prostorski povezanosti enaka: ,,Slučajno polje je prostorsko neodvisno.“ Poleg prostorsko neodvisnega slučajnega polja so lahko vrednosti v slučajnem polju urejene po neki zakonitosti. Tedaj rečemo, da so vrednosti slučajnega polja urejene. Slučajno polje je lahko tudi tako, da so vrednosti na manjšem območju bistveno drugačne od vrednosti polja na drugih delih obravnavanega območja. V tem primeru pravimo, da so vrednosti polja zbrane v gručah. Na sliki 3.13 prikazujemo primere treh polj, prostorsko neodvisnega, urejenega in polja v gručah. neodvisno urejeno v gručah SLIKA 3.13: Prostorsko neodvisno, urejeno slučajno polje in slučajno polje z vrednostmi v gručah Alternativna domneva je ,,Slučajno polje je prostorsko odvisno.“, če želimo izvesti dvostranski preizkus, ter ,,Slučajno polje je urejeno.“ ali ,,Vrednosti slučajnega polja so v gručah.“, če želimo izvesti enostranski preizkus. V geostatistiki sta se uveljavila dva preizkusa o prostorski povezanosti slučajnega polja Moranov indeks I, ki temelji na korelacijski funkciji, ter Gearyjevo razmerje C, ki temelji na variogramu. 3.8 Preizkušanje domnev o prostorski povezanosti 37 3.8.1 Moranov indeks Moranov indeks I je mera za oceno prostorske povezanosti vrednosti slučajnega polja. Pomembna razlika glede na razsevni graf, kovariančno funkcijo in variogram je, da Moranov indeks ni odvisen od razdalje med dvema točkama, temveč so v njem zajete vse razdalje med točkami. Definiran je z enačbo: n n n P P wij(Zi − ¯ Z)(Zj − ¯ Z) i=1 j=1 I = , ! (3.15) n n n P P w P ij (Zi − ¯ Z)2 i=1 j=1 i=1 kjer so z wij označene uteži, Zi in Zj so vrednosti slučajnega polja v vseh n točkah, ¯ Z pa je povprečje vrednosti slučajnega polja, ki ga izračunamo po enačbi: n ¯ 1 X Z = Zi. (3.16) n i=1 Uteži wij so lahko binarne: 1 če sta točki sosedi, wij = (3.17) 0 če točki nista sosedi. Določitev sosede je lahko vezana na razdaljo med točkami. Pogoj je torej, da je razdalja dij med sosed-njima točkama manjša od neke določene razdalje dk. Sosednjo točko lahko določimo tudi glede na skupno mejo območja, ki pripada določeni točki. Poleg tega naj velja, da točka sama sebi ni soseda, kar pomeni, da je wii = 0. Uteži so lahko tudi zvezne funkcije razdalje med točkami: wij = f (dij). (3.18) Običajno uporabimo padajoče eksponentne ali racionalne funkcije: f (d) = e−a d, 1 (3.19) f (d) = . 1 + a d Na sliki 3.14 prikazujemo dve pogosto uporabljeni funkciji, s katerima določamo vrednosti uteži. Razdaljo d običajno računamo kot Evklidovo1 razdaljo med dvema točkama, lahko pa namesto tega uporabimo tudi druge mere razdalje med točkami, kot so: razdalja Manhattan, razdalja po poteh ali cestah, čas, potreben za premik iz ene v drugo točko, stroški, potrebni za premik iz ene v drugo točko, razdalja Mahalanobis2 in drugo. 1 Evklid, grški matematik (priblližno 325–265 pred našim štetjem). 2 Prasanta Chandra Mahalanobis, indijski (bengalski) statistik (1893–1972). 38 3 Geostatistika Zvezne uteži wij f (d) f (d) = exp(-a d) f (d) = 1/(1 + a d) d SLIKA 3.14: Dve funkciji, ki ju lahko uporabimo za določitev vrednosti uteži wij Uteži lahko zapišemo v matriko uteži:  0 w  12 w13 · · · w1n  w21 0 w23 · · · w2n    W =  w31 w32 0 · · · w3n  . (3.20)  . . . .   . . . . .. .   . . . .  wn1 wn2 wn3 · · · 0 Moranov indeks I izračunamo iz geostatističnih podatkov, torej iz slučajnega vzorca, zato je slučajna spremenljivka. Zaloga vrednosti je I ∈ [−1, 1], pri čemer vrednosti okoli 0 predstavljajo prostorsko neodvisno polje, pozitiven I nakazuje na to, da so vrednosti polja zbrane v gručah, negativne vrednosti I pa nakazujejo na prostorsko pravilno urejeno polje. Srednja vrednost Moranovega indeksa je 1 E[I] = − ≈ 0, (3.21) n − 1 enačba za varianco pa je bistveno bolj zapletena: nS 2 4 − S3S5 −1 var[I] = − , (3.22) (n − 1)(n − 2)(n − 3)S2 n − 1 0 kjer Si pomenijo naslednje: n n X X S0 = wij, i=1 j=1 n n 1 X X S1 = (wij − wji)2, 2 (3.23) i=1 j=1   2 n n n X X X S2 =  wij + wji , i=1 j=1 j=1 3.8 Preizkušanje domnev o prostorski povezanosti 39 n 1 P (Z n i − ¯ Z)4 S i=1 3 = , n 2 1 P (Z n i − ¯ Z)2 i=1 (3.24) S4 = (n2 − 3n + 3)S1 − nS2 + 3S20, S5 = (n2 − n)S1 − 2nS2 + 6S20. V nadaljevanju še predpostavimo, da so geostatistični podatki taki, da je I slučajna spremenljivka I − E[I] Z = ∼ N (0, 1), (3.25) pvar[I] porazdeljena standardizirano normalno. To je osnova za preizkus domneve o prostorski neodvisnosti na osnovi Moranovega indeksa. H0: Slučajno polje je prostorsko neodvisno (I = −1/(n − 1)). H1: Slučajno polje je prostorsko odvisno (I ̸= −1/(n − 1)). α = 5 %. Zapišemo statistiko Z: I − E[I] Z = (3.26) pvar[I] in postavimo meje kritičnega območja za dvostranski preizkus domneve: Z > k1−α/2 ∨ Z < −k1−α/2 −→ H0 zavrnemo. (3.27) Če je torej Z > k1−α/2 ali Z < −k1−α/2, ničelno domnevo zavrnemo. Za α = 0.05 je kritična vrednost k0.975 = 1.96. V tem primeru izjavimo: ,,Slučajno polje je statistično značilno prostorsko odvisno!“ ali pa ,,S tveganjem, manjšim od α = 5 %, lahko trdimo, da je slučajno polje prostorsko odvisno!“ Če statistika Z ne pade v kritično območje (−k1−α/2 ≤ Z ≤ k1−α/2), potem ničelne domneve ne zavrnemo. Tedaj izjavimo: ,,Slučajno polje ni statistično značilno prostorsko odvisno!“ ali pa ,,S tveganjem, manjšim ali enakim od α = 5 %, ne moremo trditi, da je slučajno polje prostorsko odvisno!“ 40 3 Geostatistika 3.8.2 Gearyjevo razmerje Gearyjevo3 razmerje C je še ena mera za oceno prostorske povezanosti vrednosti slučajnega polja. Definirano je z naslednjo enačbo: n n P P wij(Zi − Zj)2 (n − 1) i=1 j=1 C = · , 2 ! (3.28) n n n P P w P ij (Zi − ¯ Z)2 i=1 j=1 i=1 iz katere je razvidno, da Gearyjevo razmerje izhaja iz definicije variograma. Tudi Gearyjevo razmerje je podobno kot Moranov indeks slučajna spremenljivka. Zaloga vrednosti je interval od nič do približno dve. Vrednosti C v bližini ena nakazujejo na prostorsko neodvisno slučajno polje, vrednosti C, manjše od ena, predstavljajo slučajno polje, katerega vrednosti so zbrane v gručah, medtem ko vrednosti C, večje od ena, predstavljajo prostorsko odvisna urejena polja. Srednja vrednost in varianca slučajne spremenljivke C pa sta: E[C] = 1 (3.29) in: S6 − S7 + S8 var[C] = , (3.30) n(n − 2)(n − 3)S2 0 kjer so S6 = S1(n − 1)(n2 − 3n + 3 − (n − 1)S3), 1 S7 = S2(n − 1)(n2 + 3n − 6 − (n2 − n + 2)S3), (3.31) 4 S8 = S20(n2 − 3 − (n − 1)2S3), pomen členov od S0 do S3 pa je enak kot pri Moranovemu indeksu. Podobno kot pri Moranovem indeksu predpostavimo, da so geostatistični podatki taki, da lahko tvorimo standardizirano normalno statistiko Z po enačbi: C − E[C] Z = ∼ N (0, 1). (3.32) pvar[C] To je osnova za preizkus domneve o prostorski neodvisnosti na osnovi Gearyjevega razmerja. Pokažimo sedaj primer enostranskega preizkušanja domneve o prostorski neodvisnosti. Ničelna domneva je enaka, kot pri dvostranskem preizkus, ki smo ga pokazali ob uporabi Moranovega indeksa. Pri enostranskem preizkusu v alternativni domnevi zapišemo, da so vrednosti zbrane v gručah ali pa da so vrednosti urejene po določenem vzorcu. V nadaljevanju bomo prikazali primer, ko so vrednosti zbrane v gručah. 3 Robert Charles Geary, irski statistik (1896–1983). 3.8 Preizkušanje domnev o prostorski povezanosti 41 H0: Slučajno polje je prostorsko neodvisno (C = 1). H1: Vrednosti slučajnega polja so zbrane v gručah (C < 1). α = 5 %. Zapišemo statistiko Z: C − 1 Z = (3.33) pvar[C] in postavimo meje kritičnega območja za ta enostranski preizkus domneve: Z < −k1−α −→ H0 zavrnemo. (3.34) Če je torej Z < −k1−α, ničelno domnevo zavrnemo. Za α = 0.05 je kritična vrednost k0.95 = 1.645. V tem primeru izjavimo: ,,Vrednosti so statistično značilno zbrane v gručah!“ ali pa ,,S tveganjem, manjšim od α = 5 %, lahko trdimo, da so vrednosti slučajnega polja zbrane v gručah! “ Če statistika Z ne pade v kritično območje (−k1−α ≤ Z), potem ničelne domneve ne zavrnemo. Tedaj izjavimo: ,,Vrednosti niso statistično značilno zbrane v gručah!“ ali pa ,,S tveganjem, manjšim ali enakim α = 5 %, ne moremo trditi, da so vrednosti slučajnega polja zbrane v gručah!“ Določimo vrednost Moranovega indeksa I in Gearyjevega razmerja C ter naredimo preizkus domneve o prostorski neodvisnosti za koncentracijo kroma (Cr). Nalogo bomo opravili z različno definiranimi utežmi: v primeru binarnih uteži za kritično razdaljo vzamemo dk = 1 km. V primeru zveznih uteži po enačbi (3.19) pa vzamemo, da je a = 1 km−1. Vse tri funkcije so prikazane na sliki 3.15. 1.2 f (d) = exp(-a d) 1.0 f (d) f (d) = 1/(1 + a d) 0.8 f (d) = binarno 1 ali 0 0.6 0.4 0.2 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 d SLIKA 3.15: Grafi funkcij, s katerimi smo določili uteži wij Če uporabimo binarne uteži, lahko le-te prikažemo v mreži tako, da povezavo (wij = 1) označimo s črnim kvadratkom. V enačbi (3.35) in na sliki 3.16 zapišemo matriko binarnih uteži. Na sliki prikazujemo celotno matriko za vseh 359 točk ter v povečanem merilu dva dela te matrike od točke 1 do 100 in od točke 131 do 230. 42 3 Geostatistika  0 0 1 0 0 0 0 1 0 0 · · ·  0 0 0 0 0 0 0 0 0 1 · · ·    1 0 0 0 0 0 0 1 0 0 · · ·     0 0 0 0 1 0 0 0 0 0 · · ·       0 0 0 1 0 0 0 0 0 0 · · ·    W =  0 0 0 0 0 0 0 0 0 0 · · ·  (3.35)    0 0 0 0 0 0 0 0 0 0 · · ·     1 0 1 0 0 0 0 0 0 0 · · ·     0 0 0 0 0 0 0 0 0 1 · · ·     0 1 0 0 0 0 0 0 1 0 · · ·   .  . . . . . . . . . . . . .. .. .. .. .. .. .. .. .. . . SLIKA 3.16: Matrika binarnih uteži za vseh 359 točk ter za dva dela (1–100 in 131–230) Podoben je prikaz v enačbi (3.36) in na sliki 3.17, kjer opišemo matriko uteži, ki so določene z eksponentno funkcijo (3.19). Na sliki 3.17 z različnimi stopnjami sivine prikazujemo stopnjo povezanosti med točkami. Opazimo lahko razlike v prikazu na slikah 3.17 in 3.19.  0 0.328 0.606 0.107 0.059 0.187 0.184 0.606 0.116 0.245 · · ·  0.328 0 0.247 0.171 0.130 0.071 0.102 0.205 0.299 0.430 · · ·    0.606 0.247 0 0.127 0.063 0.286 0.301 0.493 0.078 0.155 · · ·     0.107 0.171 0.127 0 0.425 0.061 0.151 0.068 0.071 0.074 · · ·       0.059 0.130 0.063 0.425 0 0.027 0.064 0.036 0.079 0.062 · · ·    W =  0.187 0.071 0.286 0.061 0.027 0 0.367 0.218 0.022 0.046 · · ·  (3.36)    0.184 0.102 0.301 0.151 0.064 0.367 0 0.160 0.031 0.053 · · ·     0.606 0.205 0.493 0.068 0.036 0.218 0.160 0 0.080 0.182 · · ·     0.116 0.299 0.078 0.071 0.079 0.022 0.031 0.080 0 0.421 · · ·     0.245 0.430 0.155 0.074 0.062 0.046 0.053 0.182 0.421 0 · · ·   .  . . . . . . . . . . . . .. .. .. .. .. .. .. .. .. . . 3.8 Preizkušanje domnev o prostorski povezanosti 43 SLIKA 3.17: Matrika uteži za eksponentno funkcijo (3.19) za vseh 359 točk ter za dva dela Za vajo lahko eno vrednost uteži, prikazano v enačbah (3.35) in (3.36), izračunamo iz podatkov v preglednici 3.1. Utež med 1. in 2. točko w12 določimo takole: koordinati 1. točke: {2.386, 3.077}, koordinati 2. točke: {2.544, 1.972}, p razdalja med točkama: d = (2.544 − 2.386)2 + (1.972 − 3.077)2 = 1.116 km, binarna utež: ker je d > 1 −→ w12 = 0, zvezna utež (eksponentna funkcija): w12 = e−1·1.116 = 0.328. Vidimo, da sta vrednosti uteži enaki kot v izrazu (3.35) oziroma (3.36). V naslednji preglednici 3.7 prikazujemo povprečje ¯ Z, vzorčno varianco S2 in statistike S0 do S8 za krom za tri različne načine določitve uteži. PREGLEDNICA 3.7: Izračun statistik za krom Statistika ¯ Z S2 S0 S1 S2 S3 Cr /binarna 35.02 113.4 22080 44160 5.75 · 106 3.316 Cr /eksponentna 35.02 113.4 26701 20551 8.22 · 106 3.316 Cr /racionalna 35.02 113.4 48340 41937 2.62 · 107 3.316 Statistika S4 S5 S6 S7 S8 Cr /binarna 5.04 · 109 4.47 · 109 2.00 · 1012 −1.53 · 1014 −1.44 · 1014 Cr /eksponentna 1.82 · 109 1.02 · 109 9.32 · 1011 −2.18 · 1014 −2.11 · 1014 Cr /racionalna 2.95 · 109 5.66 · 108 1.90 · 1012 −6.96 · 1014 −6.92 · 1014 Iz zapisanih statistik lahko izračunamo Moranov indeks I in Gearyjevo razmerje C ter oceni standardnih 44 3 Geostatistika deviacij ¯ σI in ¯ σC ter statistiki ZI in ZC. Te rezultate prikazujemo v preglednici 3.8. PREGLEDNICA 3.8: Izračun Moranovega indeksa in Gearyjevega razmerja za koncentracijo kroma Statistika I ¯ σI ZI αdej Cr /binarne 0.0732 0.00855 8.88 6.49 · 10−19 Cr /eksponentna 0.0329 0.00349 10.21 1.73 · 10−24 Cr /racionalna 0.0113 0.00146 9.62 6.86 · 10−22 Statistika C ¯ σC ZC αdej Cr /binarne 0.918 0.0214 −3.821 0.00013 Cr /eksponentna 0.952 0.0154 −3.148 0.00164 Cr /racionalna 0.979 0.00739 −2.885 0.00391 Ugotovimo lahko, da je ne glede na tip uteži wij vrednost Moranovega indeksa I mnogo višja od kritične vrednosti k1−α/2 = 1.96, zato moramo ničelno domnevo zavrniti. Dejansko tveganje, da smo sprejeli napačno odločitev, je zanemarljivo majhno (manj od 10−18). Izjavimo lahko, da je koncentracija kroma statistično značilno prostorsko odvisna. Pri preizkusu s Gearyjevim razmerjem pridemo do podobnih zaključkov, ničelno domnevo moramo zavrniti, saj je statistika ZC po absolutni vrednosti vedno večja od k1−α/2 = 1.96. Dejansko tveganje, da smo sprejeli napačno odločitev, pa je nekoliko višje kot pri preizkusu z Moranovim indeksom (od 0.0039 do 0.016), a je vseeno manjše od določenega tveganja α = 5 %. Zato tudi v tem primeru zaključimo, da je koncentracija kroma statistično značilno prostorsko odvisna. 4 Krigiranje Krigiranje je ena od pomembnejših metod, ki jih uporabljamo v prostorski statistiki. Imenuje se po južnoafriškem rudarskem inženirju D. G. Krigeju1. V različnih disciplinah se pogosto najdemo pred problemom, kjer poznamo vrednosti slučajnega polja Z(si) v končno mnogo točkah si. Naša naloga je, da ocenimo oziroma interpoliramo vrednosti polja v poljubi točki s obravnavanega območja, v kateri nimamo meritev. Interpolacijo lahko opravimo na različne načine: z linearno regresijo oziroma interpolacijo, nelinerano interpolacijo, interpolacijo z zlepki (angl. spline), interpolacijo z radialnimi baznimi funkcijami, interpolacijo z umetnimi nevronskimi mrežami in seveda, s krigiranjem. Ideja interpolacije s krigiranjem je dokaj preprosta, saj izhaja iz enačbe za uteženo povprečje vrednosti polja Z(si) v določenem številu n(s) točk v okolici točke s: n(s) Z∗(s) = m(s) + P λi(s) (Z(si) − m(si)) , (4.1) i=1 kjer je Z∗(s) interpolirana oziroma krigirana vrednost polja v točki s, λi(s) so uteži, ki jih moramo določiti za vsako interpolacijo posebej, Z(si) so znane vrednosti polja v točkah si, m(s) in m(si) so srednje vrednosti polja za točke s in si. Različni načini krigiranja se ločijo po načinu določitve srednje vrednosti m(s): preprosto, običajno krigiranje in krigiranje s trendom. Podrobnosti različnih načinov krigiranja so podane v nadaljevanju. Celoten postopek krigiranja lahko zajamemo v naslednjih točkah: 1) Iz geostatističnih podatkov Z(si) določimo vzorčni variogram (glej sliko 3.7). 2) Na osnovi vzorčnega variograma zapišemo empirično funkcijo variograma, njene koeficiente določimo po metodi najmanjših kvadratov (glej sliko 4.4). 3) Za izbrano točko s = {xK, yK} naredimo: 3.1) Poiščemo n(s) sosedov (glej sliko 4.1). 3.2) Določimo vrednosti uteži in srednjih vrednosti. 3.3) Izračunamo krigirano (interpolirano) vrednost. 3.4) Se vrnemo na 3), dokler ne interpoliramo v vseh izbranih točkah. 4) Grafično ali numerično prikažemo rezultate interpolacije. 1 Daniel Gerhardus Krige, južnoafriški rudarski inženir, avtor interpolacijske metode – krigiranja (1919–2013). 46 4 Krigiranje Na sliki 4.1 prikazujemo točko ali točke, v katerih želimo opraviti interpolacijo, ter sosede, ki so v tem primeru definirane kot prvih n(s) najbližjih točk. 6 Jura ‒ iskanje sosedov { xK,yK} = {3.3, 2.8} km 5 najbližjih 15 točk y 4 3 2 1 00 1 2 3 4 5 x SLIKA 4.1: Iskanje najbližjih točk v okolici točke, kjer želimo interpolirati vrednost polja 4.1 Preprosto krigiranje Pri preprostem krigiranju predpostavimo, da je srednja vrednost enaka na celotnem območju krigiranja: m(s) = m(si) = m. Oceno srednje vrednosti ˆ m določimo vnaprej s preprostim povprečjem vseh merjenih vrednosti slučajnega polja: n 1 X ˆ m = ¯ Z = Zi. (4.2) n i=1 Enačba krigiranja se sedaj nekoliko poenostavi: n(s) Z∗(s) = m + P λi(s) (Z(si) − m) . (4.3) i=1 Cilj postopka krigiranja je, da poiščemo take uteži λi(s), da bodo razlike med Z∗(s) in Z(s) kar najmanjše. Ta pogoj izpolnimo tako, da določimo uteži λi(s) tako, da je varianca razlike med krigirano in pravo vrednostjo najmanjša: var[Z∗(s) − Z(s)] . . . čim manjša. (4.4) 4.1 Preprosto krigiranje 47 Zato da bomo imeli pri izpeljavi enostavnejše zapise, vpeljemo novo količino R(s), ki jo imenujemo ostanek ali rezidual in predstavlja razliko med vrednostjo funkcije in njeno srednjo vrednostjo m: R(s) = Z(s) − m, R∗(s) = Z∗(s) − m, (4.5) R(si) = Z(si) − m. Iz definicije reziduala sledi, da je pričakovana vrednost reziduala enaka nič: E[R(s)] = E[Z(s) − m] = E[Z(s)] − m = m − m = 0. (4.6) Hitro lahko tudi ugotovimo, da je razlika med krigirano in pravo vrednostjo slučajnega polja enaka razliki med krigirano in pravo vrednostjo reziduala: Z∗(s) − Z(s) = (Z∗(s) − m) − (Z(s) − m) = R∗(s) − R(s), (4.7) zaradi česar je enaka tudi varianca teh razlik: var[Z∗(s) − Z(s)] = var[R∗(s) − R(s)]. (4.8) Ker so vrednosti reziduala le za konstanto spremenjene vrednosti slučajnega polja, lahko pričakujemo, da so vse lastnosti raziduala R, povezane z variabilnostjo in prostorsko odvisnostjo, enake kot za slučajno polje Z. To pomeni, da so kovariančna funkcija CR(h), variogram γR(h) in varianca CR(0) reziduala enake kovariančni funkciji C(h), variogramu γ(h) in varianci C(0) slučajnega polja: CR(h) = C(h), γR(h) = γ(h), CR(0) = C(0). (4.9) Zadnjo trditev preprosto dokažemo tako, da primerjamo kovariančni funkciji slučajnega polja in reziduala med dvema točkama si in sj. Zapišimo najprej kovariančno funkcijo slučajnega polja Z: C(si − sj) = cov[Z(si), Z(sj)] = E[(Z(si) − m)(Z(sj) − m)] = E[Z(si)Z(sj)] − m2, (4.10) podobno zapišemo še za kovariančno funkcijo reziduala R: CR(si − sj) = cov[R(si), R(sj)] = E[(R(si) − 0)(R(sj) − 0)] = E[R(si)R(sj)]. (4.11) Če v enačbo (4.11) vstavimo tretjega od izrazov (4.5), dobimo: CR(si − sj) = E[(Z(si) − m)(Z(sj) − m)] = E[Z(si)Z(sj)] − m2 = C(si − sj). (4.12) Zato bomo v nadaljevanju tudi kovariančno funkcijo reziduala označili kar s C(h). Na enak način lahko pokažemo enakost tudi za variogram in varianco. S primerjavo enačb (4.3) in (4.5) izpeljemo enačbo krigiranja za rezidual: n(s) X R∗(s) = λi(s)R(si). (4.13) i=1 48 4 Krigiranje Zapišimo varianco razlike med krigirano in pravo vrednostjo reziduala v odvisnosti od neznanih uteži λi(s) in znanih vrednosti kovariančne funkcije C(h): var[R∗(s) − R(s)] = var[R∗(s)] + var[R(s)] − 2 cov[R∗(s), R(s)], (4.14) kjer smo uporabili znano zvezo za izračun variance linearne kombinacije slučajnih spremenljivk (var[a Z1 + b Z2] = a2 var[Z1] + b2 var[Z2] + 2 ab cov[Z1, Z2]). Po vrsti bomo zapisali člene iz enačbe (4.14): var[R∗(s)] = E[R∗(s) R∗(s)] = vključimo enačbo (4.13) n(s) n(s)  X X = E  λi(s)R(si) λj(s)R(sj) = i=1 j=1 n(s) n(s)  X X = E  λi(s)λj(s)R(si)R(sj) = i=1 j=1 (4.15) upoštevamo, da je E[.] linearni operator in da uteži λi(s) niso slučajne n(s) n(s) X X = λi(s)λj(s) E [R(si)R(sj)] = i=1 j=1 zapišemo kovariančno funkcijo (enačbi (4.11) in (4.12)) in dobimo končni rezultat n(s) n(s) X X = λi(s)λj(s) C(si − sj). i=1 j=1 Z drugim členom nimamo dosti dela, saj varianco reziduala preprosto zapišemo kot kovariančno funkcijo za h = 0: var[R(s)] = C(0). (4.16) 4.1 Preprosto krigiranje 49 Tretji člen izpišemo podobno, kot prvega: 2 cov[R∗(s), R(s)] = 2 E[R∗(s) R(s)] = vključimo enačbo (4.13) n(s)  X = 2 E  λi(s)R(si)R(s) = i=1 upoštevamo, da je E[.] linearni operator in da uteži λi(s) niso slučajne (4.17) n(s) X = 2 λi(s) E [R(si)R(s)] = i=1 zapišemo kovariančno funkcijo (enačbi (4.11) in (4.12)) in dobimo končni rezultat n(s) X = 2 λi(s) C(si − s). i=1 Varianca razlike krigiranega in pravega reziduala je odvisna le od neznanih uteži λi(s) ter znanih vrednosti kovariančne funkcije C(si − sj) in C(si − s) in variance C(0): n(s) n(s) n(s) X X X var[R∗(s) − R(s)] = λi(s)λj(s) C(si − sj) + C(0) − 2 λi(s) C(si − s). (4.18) i=1 j=1 i=1 Z odvajanjem enačbe (4.18) po utežeh λi(s) lahko določimo take vrednosti uteži λi(s), da bo varianca najmanjša: n(s) ∂(var[R∗(s) − R(s)]) X = 2 λj(s)C(sk − sj) − 2 C(sk − s) = 0, (4.19) ∂λk(s) j=1 kjer je k = 1, . . . , n(s). Zgornja enačba predstavlja sistem n(s) linearnih enačb za prav toliko neznank λi(s). Izpišimo nekaj enačb tega sistema v razviti obliki: C(0)λ1(s) + C(s1 − s2)λ2(s) + · · · + C(s1 − sn(s))λn(s)(s) = C(s1 − s), C(s2 − s1)λ1(s) + C(0)λ2(s) + · · · + C(s2 − sn(s))λn(s)(s) = C(s2 − s), .. . . . . . .. . . .. .. C(sn(s) − s1)λ1(s) + C(sn(s) − s2)λ2(s) + · · · + C(0)λn(s)(s) = C(sn(s) − s). Sistem linearnih enačb običajno zapišemo v matrični obliki:  C(0) C(s      1 − s2) · · · C(s1 − sn(s)) λ1(s) C(s1 − s)  C(s2 − s1) C(0) · · · C(s2 − sn(s))   λ2(s)   C(s2 − s)   . . .   .  =  .  .  . . . .. .   ..   ..   . . .      C(sn(s) − s1) C(sn(s) − s2) · · · C(0) λn(s)(s) C(sn(s) − s) (4.20) 50 4 Krigiranje Z rešitvijo sistema enačb (4.20) določimo vrednosti uteži, ki jih uporabimo pri izračunu interpolacije slučajnega polja po enačbi (4.3). 4.2 Običajno krigiranje Pri običajnem krigiranju predpostavimo, da je v okolici točke, za katero interpoliramo vrednost, srednja vrednost m(s) konstantna, a neznana. Zato predpostavimo: m(s) = m(si) = µ(s). (4.21) Enačbo krigiranja lahko zapišemo takole: n(s) X Z∗(s) = µ(s) + λi(s) (Z(si) − µ(s)) , i=1 po manjši preureditvi pa v naslednji obliki: n(s) ! n(s) Z∗(s) = µ(s) 1 − P λ P i(s) + λi(s)Z(si). (4.22) i=1 i=1 Ker ne želimo, da bi izbira srednje vrednosti vplivala na krigirano vrednost, zahtevamo, da je prvi člen v enačbi (4.22) enak nič. To lahko dosežemo le, če je izpolnjen pogoj: n(s) X λi(s) = 1. (4.23) i=1 Prišli smo do problema vezanih ekstremov, kjer mora biti ob določitvi minimuma var[Z∗(s) − Z(s)] izpolnjen tudi pogoj (4.23). Problem rešimo z uvedbo Lagrangevega2 množitelja 2µ(s): n(s)  X L(λi(s), µ(s)) = var[Z∗(s) − Z(s)] + 2µ(s)  λi(s) − 1 (4.24) i=1 oziroma v razviti obliki: n(s) n(s) X X L(λi(s), µ(s)) = λi(s)λj(s) C(si − sj) + C(0)− i=1 j=1 (4.25) n(s) n(s)  X X − 2 λi(s) C(si − s) + 2µ(s)  λi(s) − 1 . i=1 i=1 2 Joseph-Louis Lagrange, francoski matematik, rojen v Torinu (1736–1813). 4.3 Krigiranje s trendom 51 Z odvajanjem L(λi(s), µ(s)) po λi(s) in µ(s) dobimo sistem n(s) + 1 lineranih enačb za prav toliko neznank: ∂L(λi(s), µ(s)) = 0, ∂λk(s) (4.26) ∂L(λi(s), µ(s)) = 0. ∂µ(s) Izpišimo nekaj značilnih enačb sistema, pri čemer zaradi preglednejšega zapisa upoštevamo λi(s) = λi in µ(s) = µ: C(0) λ1 + C(s1 − s2) λ2 + · · · + C(s1 − sn(s)) λn(s) + 1 µ = C(s1 − s), C(s2 − s1) λ1 + C(0) λ2 + · · · + C(s2 − sn(s)) λn(s) + 1 µ = C(s2 − s), .. . . . . . .. . . .. .. C(sn(s) − s1) λ1 + C(sn(s) − s2) λ2 + · · · + C(0) λn(s) + 1 µ = C(sn(s) − s), 1 µ λ1 + 1 λ2 + · · · + 1 λn(s) + 0 µ = 1. Sistem linearnih enačb običajno zapišemo v matrični obliki:  C(0) C(s      1 − s2) · · · C(s1 − sn(s)) 1 λ1 C(s1 − s)  C(s2 − s1) C(0) · · · C(s2 − sn(s)) 1   λ2   C(s2 − s)   . . . . .   .   .   . . . . . .   .  =  .  .  . . . .   .   .     λ   C(s   C(sn(s) − s1) C(sn(s) − s2) · · · C(0) 1   n(s)   n(s) − s)  1 1 · · · 1 0 µ 1 (4.27) Ker z dodatno enačbo izpolnimo pogoj (4.23), pri določitvi krigirane vrednosti Z∗(s) po enačbi (4.22) lahko izpustimo prvi člen in se enačba krigiranje še poenostavi: n(s) Z∗(s) = P λi(s)Z(si). (4.28) i=1 4.3 Krigiranje s trendom Pri krigiranju s trendom predpostavimo, da je v okolici točke, za katero interpoliramo vrednost, srednja vrednost m(s) oziroma m(si) znana funkcija z neznanimi koeficienti. Zato predpostavimo: K K X X m(s) = αk(s)fk(s) oziroma m(si) = αk(s)fk(si), (4.29) k=1 k=1 kjer je K število neznanih koeficientov αk(s) oziroma število funkcij fk(s), ki sestavljajo funkcijo trenda. Koeficiente trenda določimo v postopku krigiranja skupaj z utežmi λi(s). Enačbo krigiranja 52 4 Krigiranje lahko zapišemo takole: K n(s) K ! X X X Z∗(s) = αk(s)fk(s) + λi(s) Z(si) − αk(s)fk(si) , k=1 i=1 k=1 po manjši preureditvi pa v naslednji obliki:   K n(s) n(s) X X X Z∗(s) = αk(s) fk(s) − λi(s)fk(si)+ λi(s)Z(si). (4.30) k=1 i=1 i=1 Ker ne želimo, da bi izbira funkcije trenda vplivala na krigirano vrednost, zahtevamo, da je prvi člen v enačbi (4.30) enak nič. To lahko dosežemo le, če so izpolnjeni pogoji: n(s) X λi(s)fk(si) = fk(s), (4.31) i=1 za k = 1, . . . , K. Prišli smo do problema vezanih ekstremov, kjer morajo biti ob določitvi minimuma var[Z∗(s) − Z(s)] izpolnjeni tudi pogoji (4.31). Problem rešimo z uvedbo Lagrangevih množiteljev 2αk(s):   K n(s) X X L(λi(s), αk(s)) = var[Z∗(s) − Z(s)] + 2 αk(s)  λi(s)fk(si) − fk(s) (4.32) k=1 i=1 oziroma v razviti obliki: n(s) n(s) X X L(λi(s), αk(s)) = λi(s)λj(s) C(si − sj) + C(0)− i=1 j=1 n(s) X − 2 λi(s) C(si − s)+ (4.33) i=1   K n(s) X X + 2 αk(s) λ .  i(s)fk(si) − fk(s) k=1 i=1 Z odvajanjem L(λi(s), αk(s)) po λi(s) in αk(s) dobimo sistem n(s) + K lineranih enačb za prav toliko neznank: ∂L(λi(s), αk(s)) = 0, ∂λl(s) (4.34) ∂L(λi(s), αk(s)) = 0. ∂αm(s) 4.3 Krigiranje s trendom 53 Izberemo lahko različne enačbe trenda. Pogosto uporabimo najpreprostejši linearni trend, ki ga za dve dimenziji zapišemo z naslednjo enačbo: m(s) = m(x, y) = a1 + a2 x + a3 y. (4.35) V primeru, da želimo uporabiti kvadratni trend, je enačba trenda enaka: m(s) = m(x, y) = a1 + a2 x + a3 y + a4 x2 + a5 y2 + a6 x y. (4.36) Pri linearnem trendu v dveh dimenzijah moramo določiti tri neznane koeficiente K = 3, pri kvadratnem trendu pa je K = 6. Na sliki 4.2 prikazujemo tri različne funkcije trenda, ki bi jih lahko uporabili pri krigiranju s trendom: linearni, kvadratni in kubični trend. V primeru, da izberemo konstantno funkcijo trenda: m(s) = m(x, y) = a1, (4.37) se primer poenostavi v običajno krigiranje, ki smo ga obravnavali v prejšnjem poglavju. SLIKA 4.2: Primeri funkcij trenda (linearna, kvadratna in kubična funkcija) Izpišimo nekaj značilnih enačb sistema za linearno enačbo trenda (4.35), pri čemer zaradi preglednejšega zapisa upoštevamo λi(s) = λi in αk(s) = αk: C(0) λ1 + · · · + C(s1 − sn(s)) λn(s) + 1 α1 + x1 α2 + y1 α3 = C(s1 − s), C(s2 − s1) λ1 + · · · + C(s2 − sn(s)) λn(s) + 1 α1 + x2 α2 + y2 α3 = C(s2 − s), .. . . . . . . . . . .. .. .. .. .. C(sn(s) − s1) λ1 + · · · + C(0) λn(s) + 1 α1 + xn(s) α2 + yn(s) α3 = C(sn(s) − s), 1 λ1 + · · · + 1 λn(s) + 0 α1 + 0 α2 + 0 α3 = 1, x1 λ1 + · · · + xn(s) λn(s) + 0 α1 + 0 α2 + 0 α3 = x, y1 λ1 + · · · + yn(s) λn(s) + 0 α1 + 0 α2 + 0 α3 = y. 54 4 Krigiranje Sistem linearnih enačb običajno zapišemo v matrični obliki:  C(0) C(s      1 − s2) · · · C(s1 − sn(s)) 1 x1 y1 λ1 C(s1 − s) C(s λ C(s  2 − s1) C(0) · · · C(s2 − sn(s)) 1 x2 y2   2   2 − s)   . . . . . . .   .   .   .. .. . . .. .. .. ..   ..   ..         C(s   λ  =  C(s  .  n(s) − s1) C(sn(s) − s2) · · · C(0) 1 xn(s) yn(s)   n(s)   n(s) − s)         1 1 · · · 1 0 0 0   α1   1         x1 x2 · · · xn(s) 0 0 0   α2   x  y1 y2 · · · yn(s) 0 0 0 α3 y (4.38) Ker z dodatnimi K enačbami zagotovimo, da so izpolnjeni pogoji (4.31), lahko prvi člen v enačbi (4.30) izpustimo in se enačba krigiranja še poenostavi: n(s) Z∗(s) = P λi(s)Z(si). (4.39) i=1 S primerjavo enačb (4.20), (4.27) in (4.38) vidimo, da prvih n(s) enačb v (4.27) in (4.38) predstavlja sistem preprostega krigiranja, prvih n(s) + 1 enačb v (4.38) pa predstavlja sistem enačb običajnega krigiranja. Različni deli sistema so v enačbi (4.38) ločeni z rdečimi črtami. 4.4 Primeri krigiranja Vrnimo se na že večkrat obravnavani primer geostatističnih podatkov z območja Jure. Ocenimo vrednost koncentracije kroma v točki Tk (xk = 3.3 km in yk = 2.8 km), ki leži približno na sredini obravnavanega območja. Točka s pripadajočimi točkami v okolici je prikazana na sliki 4.1. Na podlagi podatkov o koncentraciji kroma smo ocenili variogram, pripravili empirični variogram in ustrezno kovariančno funkcijo, ki sta prikazana na sliki 4.3. Pri aproksimaciji variograma smo uporabili enačbo (3.7), vrednosti koeficientov pa so prikazane v preglednici 3.5. Pri določitvi empirične kovariančne funkcije, ki jo potrebujemo pri krigiranju, smo uporabili zvezo: C(h) = C(0) − γ(h) = 77.3 e−4.97h. (4.40) 4.4 Primeri krigiranja 55 120 120 C(h) g(h) 80 80 40 40 0 0 0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 h h SLIKA 4.3: Empirični variogram in ustrezna kovariančna funkcija V krigiranje vključimo deset najbližjih točk v okolici točke Tk. Koordinate teh točk in koncentracije kroma v teh točkah so prikazane v preglednici 4.1. Za izračun uteži λi moramo izračunati razdalje med temi desetimi točkami, razdalje med temi desetimi točkami in točko Tk ter uporabiti enačbo empirične kovariančne funkcije (4.40). 4.4.1 Preprosto krigiranje Pri preprostem krigiranju moramo oceno srednje vrednosti določiti vnaprej: n 1 X m = ¯ Z = Zi = 35.02. (4.41) n i=1 PREGLEDNICA 4.1: Koordinate in vrednosti koncentracije kroma v 10 najbližjih točkah točki Tk xi yi Zi 3.316 2.987 33.00 3.077 2.926 35.32 3.312 2.991 32.60 3.347 2.505 33.20 3.322 3.001 35.00 3.308 2.963 36.00 3.422 2.851 38.72 3.557 2.640 25.72 3.212 2.716 35.52 3.287 3.061 45.92 56 4 Krigiranje Sistem enačb za določitev uteži pri preprostem krigiranju je:  77.3 22.7 75.2 7.0 71.7 68.2 32.8 9.5 18.3 52.1   λ    1 30.4 22.7 77.3 23.0 6.4 21.6 24.2 13.4 4.8 22.4 22.4 λ 21.7    2     75.2 23.0 77.3 6.9 72.1 67.2 31.9 9.2 18.1 53.4   λ   29.9     3     7.0 6.4 6.9 77.3 6.6 7.9 13.3 22.4 22.3 4.8   λ   17.5     4           71.7 21.6 72.1 6.6 77.3 63.2 31.6 9.1 16.9 54.7   λ5   28.3      =   ,  68.2 24.2 67.2 7.9 63.2 77.3 34.9 10.2 20.7 47.0   λ6   34.4         32.8 13.4 31.9 13.3 31.6 34.9 77.3 22.3 22.4 22.4   λ7   40.1         9.5 4.8 9.2 22.4 9.1 10.2 22.3 77.3 13.4 6.4   λ8   17.2         18.3 22.4 18.1 22.3 16.9 20.7 22.4 13.4 77.3 13.4   λ9   42.2  52.1 22.4 53.4 4.8 54.7 47.0 22.4 6.4 13.4 77.3 λ10 21.1 rešitev tega sistema pa je:  λ    1 −0.0102 λ 0.0586  2     λ   0.0026   3     λ   0.0327   4       −   λ5   0.0224    =   .  λ6   0.2103       λ7   0.3036       λ8   0.0328       λ9   0.3784  λ10 −0.0090 Vrednost koncentracije kroma v točki Tk pri preprostem krigiranju izračunamo po enačbi (4.3): Z∗(3.3, 2.8) = 35.02 + (−0.0102 (33.00 − 35.02) + 0.0586 (35.32 − 35.02) + · · · ) = 36.11. 4.4.2 Običajno krigiranje Na podoben način izračunamo tudi uteži in oceno srednje vrednosti pri običajnem krigiranju. Sistem enačb se poveča za eno enačbo in eno neznanko µ:       77.3 22.7 75.2 7.0 71.7 68.2 32.8 9.5 18.3 52.1 1.0 λ1 30.4  22.7 77.3 23.0 6.4 21.6 24.2 13.4 4.8 22.4 22.4 1.0   λ2   21.7         75.2 23.0 77.3 6.9 72.1 67.2 31.9 9.2 18.1 53.4 1.0   λ3   29.9         7.0 6.4 6.9 77.3 6.6 7.9 13.3 22.4 22.3 4.8 1.0   λ4   17.5         71.7 21.6 72.1 6.6 77.3 63.2 31.6 9.1 16.9 54.7 1.0   λ5   28.3         68.2 24.2 67.2 7.9 63.2 77.3 34.9 10.2 20.7 47.0 1.0   λ    6 = 34.4 ,        32.8 13.4 31.9 13.3 31.6 34.9 77.3 22.3 22.4 22.4 1.0   λ    7 40.1        9.5 4.8 9.2 22.4 9.1 10.2 22.3 77.3 13.4 6.4 1.0   λ   17.2     8     18.3 22.4 18.1 22.3 16.9 20.7 22.4 13.4 77.3 13.4 1.0   λ   42.2     9     52.1 22.4 53.4 4.8 54.7 47.0 22.4 6.4 13.4 77.3 1.0   λ   21.1     10    1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 µ 1.0 4.4 Primeri krigiranja 57 rešitev tega sistema pa je:     λ1 −0.0099  λ2   0.0628       λ3   0.0021       λ4   0.0370       λ5   −0.0211       λ    6 = 0.2106 .      λ    7 0.3060      λ   0.0371   8     λ   0.3809   9     λ   −0.0055   10    µ −0.5671 Vrednost koncentracije kroma v točki Tk pri običajnem krigiranju izračunamo po enačbi (4.28): Z∗(3.3, 2.8) = (−0.0099 · 33.00 + 0.0628 · 35.32 + 0.0021 · 32.60 · · · ) = 36.11. Vidimo, da smo prišli do enakega rezultata kot pri preprostem krigiranju, kar pa ni pravilo. Pri različnih točkah oziroma podatkih se rezultati preprostega in običajnega krigiranja lahko razlikujejo. 4.4.3 Krigiranje s trendom Ponovimo izračun še za krigiranje s trendom, pri čemer se odločimo za linearni trend. Sistem enačb se še nekoliko razˇsiri:3 77.3 22.7 75.2 7.0 71.7 68.2 32.8 9.5 18.3 52.1 1.0 3.3 3.0  λ    1 30.4 22.7 77.3 23.0 6.4 21.6 24.2 13.4 4.8 22.4 22.4 1.0 3.1 2.9  λ2  21.7       75.2 23.0 77.3 6.9 72.1 67.2 31.9 9.2 18.1 53.4 1.0 3.3 3.0  λ3  29.9        7.0 6.4 6.9 77.3 6.6 7.9 13.3 22.4 22.3 4.8 1.0 3.3 2.5  λ4  17.5       71.7 21.6 72.1 6.6 77.3 63.2 31.6 9.1 16.9 54.7 1.0 3.3 3.0  λ5  28.3       68.2 24.2 67.2 7.9 63.2 77.3 34.9 10.2 20.7 47.0 1.0 3.3 3.0  λ6  34.4       32.8 13.4 31.9 13.3 31.6 34.9 77.3 22.3 22.4 22.4 1.0 3.4 2.9  λ    7 = 40.1 ,        9.5 4.8 9.2 22.4 9.1 10.2 22.3 77.3 13.4 6.4 1.0 3.6 2.6  λ    8 17.2       18.3 22.4 18.1 22.3 16.9 20.7 22.4 13.4 77.3 13.4 1.0 3.2 2.7  λ  42.2    9    52.1 22.4 53.4 4.8 54.7 47.0 22.4 6.4 13.4 77.3 1.0 3.3 3.1 λ  21.1    10    1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0  α   1.0     1     3.3 3.1 3.3 3.3 3.3 3.3 3.4 3.6 3.2 3.3 0.0 0.0 0.0  α   3.3     2    3.0 2.9 3.0 2.5 3.0 3.0 2.9 2.6 2.7 3.1 0.0 0.0 0.0 α3 2.8 3Elementi matrike so zaradi omejenega prostora zapisane le na eno decimalko natančno, zaradi česar se izgubijo prave vrednosti koordinat točk. V izračunu smo upoštevali natančnejše vrednosti. 58 4 Krigiranje rešitev tega sistema pa je:  λ    1 −0.0105  λ2   0.0665       λ3   0.0032       λ4   0.0414       λ5   −0.0233       λ6   0.2108       λ    7 = 0.3030 .      λ    8 0.0335      λ   0.3843   9     λ   −0.0088   10     α   −9.4340   1     α   1.6570   2    α3 1.2100 Vrednost koncentracije kroma v točki Tk pri krigiranju s trendom izračunamo po enačbi (4.39): Z∗(3.3, 2.8) = (−0.0105 · 33.00 + 0.0665 · 35.32 + 0.0032 · 32.60 · · · ) = 36.09, ki je le malo drugačna od vrednosti, ki smo jo pridobili s preprostim oziroma običajnim krigiranjem. Pri modeliranju neznanih vrednosti slučajnega polja najpogosteje določamo aproksimacijo vzdolže neke osi ali po nekem območju. Na sliki 4.4 prikazujemo rezultate krigiranja vzdolž osi y pri x = 3.15 km. Uporabili smo vse tri vrste krigiranja: preprosto (modra), običajno (rdeča) in krigiranje z linearnim trendom (črtkana zelena). Poleg tega smo krigirali z različnim številom okoliških točk, od n(s) = 6 do n(s) = 42. Opazimo lahko, da so razlike med preprostim in običajnim krigiranjem zelo majhne, pa še to le pri manjšem številu okoliških točk. Krigiranje s trendom se nekoliko bolj razlikuje še posebej za manjše število okoliških točk in na obeh robovih območja, kjer je lahko vpliv trenda večji. 4.4 Primeri krigiranja 59 6 Izbira točk vzdolž osi y n(s) = 15 5 4 100 Krigiranje vzdolž osi y y 3 80 n(s) = 6 60 2 Cr 40 1 20 00 1 2 3 4 5 6 00 1 2 3 4 5 y x SLIKA 4.4: Krigiranje vzdolž osi y pri x = 3.15 km Na sliki 4.5 prikazujemo rezultate krigiranja za celotno območje Jura pri preprostem in običajnem krigiranju za n(s) = 10. Vidimo, da krigiranje sledi spreminjanju vrednosti slučajnega polja pri obeh načinih krigiranja. Na spodnji sliki vidimo, da se rezultati preprostega in običajnega krigiranja na nekaterih mestih razlikujejo. preprosto krigiranje običajno krigiranje y y x x y x primerjava SLIKA 4.5: Rezultati krigiranja koncentracije kroma na celotnem območju Jura 60 4 Krigiranje 5 Prostorsko utežena regresija S pojmom regresija opisujemo postopek za določitev ocen parametrov določenih matematičnih modelov. Če pojav Z opisujemo z linearnim modelom f , ki ga v k-dimenzionalnem prostoru zapišemo z naslednjo enačbo: z∗(x1, x2, . . . , xk) = f (x1, x2, . . . , xk) = a0 + a1 x1 + a2 x2 + · · · + ak xk, (5.1) potem z regresijsko analizo določimo ocene parametrov âi, i = 0, 1, 2 . . . k. Taki regresiji pravimo linearna regresija. Podatke, ki jih potrebujemo za regresijo, lahko zapišemo v naslednji preglednici in imajo enako strukturo kot geostatistični podatki (glej preglednico 3.1). PREGLEDNICA 5.1: Struktura podatkov za regresijo X1 X2 · · · Xk Z x11 x12 · · · x1k z1 x21 x22 · · · x2k z2 .. . . . . . .. . . .. .. xn1 xn2 · · · xnk zn 5.1 Linearna regresija Opišimo najprej preprost primer linearne regresije za eno samo neodvisno spremenljivko x. V tem primeru se enačba (5.1) poenostavi: z∗(x) = f (x) = a0 + a1 x. (5.2) Na sliki 5.1 je prikazan primer linearne regresije za eno neodvisno spremenljivko. Označena sta tudi ena od točk (xi, yi) in odstopanje oziroma rezidual εi. Ordinata zi je vsota vrednosti funkcije f (xi) in reziduala εi: zi = f (xi) + εi −→ εi = zi − f (xi) = zi − (a0 + a1 xi). (5.3) Za določitev ocen â0 in â1 parametrov a0 in a1 postavimo pogoj, da bi bila odstopanja εi čim manjša. Ob predpostavki, da se odstopanja porazdeljujejo normalno s srednjo vrednostjo nič, je najprimernejša 62 5 Prostorsko utežena regresija metoda najmanjših kvadratov, pri kateri postavimo pogoj, da je vsota kvadratov odstopanj čim manjša: n n n X X X S = ε2i = (zi − f (xi))2 = (zi − â0 − â1 xi)2 naj bo čim manjša. (5.4) i=1 i=1 i=1 6 5 z 4 xi 3 εi 2 z f (x i i) 1 00 1 2 3 4 5 6 x SLIKA 5.1: Linearna regresija za eno neodvisno spremenljivko V primeru, da se odstopanja ne porazdeljujejo normalno oziroma vključujejo osamelce (angl. outlier), pa je bolje uporabiti eno izmed različic robustne regresijske metode: absolutne vrednosti: n X S = |εi| naj bo čim manjša, (5.5) i=1 Huberjeva1 metoda: n X 1/2 ε2 |ε S = i i| ≤ h naj bo čim manjša, (5.6) h|εi| − 1/2 h2 |εi| > h i=1 bikvadratna enačba: n X h2/6 (1 − (1 − (ε S = i/h)2)3) |εi| ≤ h naj bo čim manjša, (5.7) h2/6 |εi| > h i=1 kjer s h označimo velikost območja, na katerem je funkcija S nelinearno odvisna od εi. V nadaljevanju se bomo omejili le na metodo najmanjših kvadratov. Iz sistema linearnih enačb, ki jih izpeljemo iz pogoja (5.4):  n   n  X   X n x â i 0 zi      i=1     i=1   n n    =  n  , (5.8)  X X     X   xi x2i  â  z  1 i xi i=1 i=1 i=1 1 Peter Jost Huber, švicarski statistik, utemeljitelj robustne statistike (1934–). 5.1 Linearna regresija 63 določimo vrednosti ocen â0 in â1: SXZ SXZ â ¯ 0 = ¯ Z − X, â1 = , (5.9) S2 S2 X X kjer sta ¯ X in ¯ Z povprečji vzorcev xi in zi, SXZ je vzorčna kovarianca med X in Z, S2 pa je vzorčna X varianca X. V primeru, da izvajamo linearno regresijo za k spremenljivk, iz sistema normalnih enačb  n P x P     P  i1 xi2 · · · P xik â0 zi P P P P  xi1 x2 x â z i1 i2 xi1 · · · P xik xi1   1   i xi1   P P P     P   xi2 xi1 xi2 x2 · · · P x â z i2 ik xi2   2  =  i xi2  (5.10)  . . . .   .   .   . . . . .. .   .   .   . . . .   .   .  P x P P P ik xi1 xik xi2 xik · · · P x2 â z ik k i xik izračunamo neznane ocene parametrov â0, â1, . . . , âk, ki so enake na celotnem obravnavanem območju. Linearno regresijo več spremenljivk lahko zapišemo tudi v matričnem zapisu. Definiramo naslednje matrike: matriko podatkov X, matriko neznanih koeficientov a, matriko vrednosti odvisnih spremenljivk Y in matriko odstopanj e. Matrike zapišemo z naslednjimi izrazi:  1 x        11 x12 · · · x1k a0 z1 ε1  1 x21 x22 · · · x2k   a1   z2   ε2          X =  1 x31 x32 · · · x3k  , a =  a2  , Z =  z3  , e =  ε3  .  . . . .   .   .   .   . . . . .. .   .   .   .   . . . .   .   .   .  1 xn1 xn2 · · · xnk ak zn εn Linearni model zapišemo z enačbo: Z = Xa + e, (5.11) iz katere lahko izpeljemo enačbo odstopanj: e = Z − Xa. Po metodi najmanjših kvadratov zapišemo najprej vsoto kvadratov odstopanj S: S = eT e = (Z − Xa)T (Z − Xa) . (5.12) Minimum vsote kvadratov odstopanj dobimo z odvajanjem S po neznanih koeficientih a in upoštevaje, da je ta odvod enak nič: dS = 0 → 2 XT (Z − Xâ) = 2 XT Z − XT Xâ = 0, da (5.13) → XT Xâ = XT Z, kar predstavlja sistem normalnih enačb v matričnem zapisu. 64 5 Prostorsko utežena regresija 5.2 Utežena linearna regresija Zaradi prostorske odvisnosti, nestacionarnosti in drugih vplivov parametri a0, a1, . . . niso enaki na celotnem obravnavanem območju. To je osnovna ideja utežene linearne regresije, pri kateri za vsako točko, v kateri želimo določiti regresijo, izračunamo drugačne ocene parametrov. Razlika izhaja iz spremenjene enačbe metode najmanjših kvadratov (5.4), ki ji dodamo uteži wi: n X S = wi ε2i naj bo čim manjša. (5.14) i=1 Ob upoštevanju, da uteži wi uredimo v matriko W:  w  1 0 · · · 0  0 w2 · · · 0  W =  . . . .  , (5.15)  .. .. . . ..    0 0 · · · wn vsoto kvadratov odstopanj zapišemo takole: S = eT W e = (Z − Xa)T W (Z − Xa) . (5.16) Minimum vsote kvadratov odstopanj dobimo z odvajanjem S po neznanih koeficientih a in zahtevo, da je ta odvod enak nič: dS = 0 → 2 XT W (Z − Xâ) = 2 XT W Z − XT W Xâ = 0, da (5.17) → XT W X â = XT W Z, kar predstavlja sistem normalnih enačb v matričnem zapisu ob upoštevanju različno uteženih odstopanj. 5.2.1 Določitev uteži Poseben izziv pri uteženi linearni regresiji je določitev uteži wi. Intuitivno bomo predpostavili, da naj vrednosti v točkah, ki so bližje točki T s koordinatami (xT 1, xT 2, . . . , xT k), bolj vplivajo na izračunano vrednost z∗ kot vrednosti v bolj oddaljenih točkah. Uteži wi določimo s padajočo funkcijo razdalje di med točkami Ti, kjer poznamo vrednosti Zi, in točko T , za katero vrednost določimo z modelom (5.1). Razdaljo di običajno izračunamo kot Evklidovo razdaljo med točkama T in Ti: q p di = (xT 1 − xi1)2 + (xT 2 − xi2)2 + · · · + (xT k − xik)2 = (XT − Xi)T (XT − Xi). (5.18) Vrednosti uteži lahko določimo z eno od naslednjih enačb: 1 d ≤ h w(d) = , (5.19) 0 d > h 5.2 Utežena linearna regresija 65 ( 1 − d 22 d ≤ h w(d) = h , (5.20) 0 d > h w(d) = e−( d )2 h . (5.21) Enačba (5.19) predstavlja binarno odvisnost, kjer je vrednost uteži enaka 1, če je razdalja manjša od kritične razdalje oziroma praga h. Pri enačbi (5.20) vrednost uteži w zvezno pada proti nič, ki jo doseže pri kritični razdalji d = h. Pri enačbi (5.21) vrednost uteži w zvezno pada proti nič, a je nikoli ne doseže. Kritična razdalja h pri tej enačbi predstavlja parameter, ki pove, kako hitro pada vrednost uteži. 1.2 enačba (6.19) 1.0 w enačba (6.20) i 0.8 enačba (6.21) 0.6 0.4 0.2 0.0 0 1 2 3 4 di SLIKA 5.2: Vrednosti uteži v odvisnosti od razdalje med točkami 5.2.2 Izbira funkcije za določitev uteži in kritične razdalje V prejšnjem razdelku smo pojasnili, kako določimo uteži v odvisnosti od razdalje točke, kjer poznamo vrednost slučajnega polja, in točke, za katero želimo to vrednost oceniti. Prikazali smo tri različne funkcije (5.19)–(5.21), ki so odvisne tudi od parametra h. Katera funkcija in kolikšen parameter h je optimalna rešitev, moramo še določiti. Ena od možnosti je izračun vrednosti navzkrižne validacije (angl. crossvalidation score), pri kateri opravimo prostorsko uteženo regresijo za različne funkcije in različne vrednosti h v vseh točkah, kjer poznamo vrednost slučajnega polja. Pri vsakokratni določitvi utežene regresije s seznama točk s poznanimi vrednostmi slučajnega polja izključimo prav točko, za katero določamo uteženo regresijo. V tem primeru torej določamo regresijo z n−1 točk, kjer je n število točk, v katerih poznamo vrednosti slučajnega polja. Tako določene ocene vrednosti slučajnega polja primerjamo z izmerjenimi oziroma tvorimo naslednjo cenilko: n X CV (h) = (zi − z∗i(h))2. (5.22) i=1 Najnižja vrednost CV (h) ustreza optimalni funkciji in vrednosti h. Izračun vrednosti CV (h) za različne enačbe bomo prikazali v naslednjem razdelku. 66 5 Prostorsko utežena regresija 5.2.3 Primeri Vzemimo že znane podatke z območja Jura. Za vrednosti koncentracije kroma opravimo običajno linearno regresijo in uteženo regresijo za vse tri tipe utežnih funkcij (enačbe (5.19)–(5.21)) pri določenih kritičnih razdaljah h. Model linearne regresije za celotno območje je: Z(x, y) = 37.0827 − 0.842091x + 0.163621y. (5.23) V preglednici 5.2 prikazujemo izračun razdalj med desetimi najbližjimi točkami, kjer poznamo vrednosti, in točko T , ki je podana s koordinatama xT = 2 km in yT = 2 km. V isti preglednici so tudi uteži, izračunane po enačbah (5.19)–(5.21), s kritičnimi razdaljami h = 0.3 km za prvi dve enačbi in h = 0.1 km za zadnjo. PREGLEDNICA 5.2: Izračun uteži za nekaj točk za točko T = (2, 2) km x y Z d w (5.19) w (5.20) w (5.21) [km] [km] / [km] / / / 1.948 1.906 49.60 0.107 1 0.7600 0.3154 2.159 2.041 37.72 0.164 1 0.4906 0.06746 1.813 2.116 33.48 0.220 1 0.2134 0.007887 2.024 2.251 38.00 0.252 1 0.0862 0.001733 2.083 1.695 29.28 0.316 0 0 4.579·10−5 1.720 1.833 31.68 0.326 0 0 2.421·10−5 1.709 1.843 33.32 0.331 0 0 1.786·10−5 1.706 1.848 34.84 0.331 0 0 1.749·10−5 2.294 1.830 46.80 0.340 0 0 9.795·10−6 1.690 1.858 38.72 0.341 0 0 8.927·10−6 Ob upoštevanju vseh 359 točk, kjer so podane vrednosti slučajnega polja, lahko po enačbi (5.17) ocenimo parametra a0 in a1 ter izračunamo vrednost polja v točki T. V preglednici 5.3 prikazujemo rezultate običajne in utežene regresije v točki T = (2, 2) km. PREGLEDNICA 5.3: Vrednost Z∗ v točki T = (2, 2) km z različnimi linearnimi regresijami Običajna Utežena Utežena Utežena regresija regresija regresija regresija (5.19) (5.20) (5.21) 35.73 42.49 42.80 42.52 Vidimo, da je rezultat utežene regresije v vseh treh primerih višji od vrednosti, ki smo jo dobili z običajno regresijo. To je rezultat dejstva, da so vrednosti v okolici točke T = (2, 2) km izmerjene vrednosti 5.2 Utežena linearna regresija 67 slučajnega polja višje od vrednosti, ki jo dobimo s preprosto linearno regresijo, kar je razvidno tudi iz preglednice 5.2. Na sliki 5.3 prikazujemo tri regresijske ploskve, ki smo jih določili z uteženo linearno regresijo (5.17) s tremi različnimi enačbami za uteži, z rumeno ploskvijo pa je prikazan rezultat običajne linearne regresije. Na spodnji desni sliki na istem grafu primerjamo vse tri različice utežene linearne regresije in običajno linearno regresijo. Opazimo lahko, da različni načini določitve uteži pripeljejo do različnih ploskev, ki bolj ali manj intenzivno sledijo podatkom, ki odstopajo od linearne ploskve, ki jo dobimo z običajno linearno regresijo. Uteži po enačbi (5.19) h = 0.3 Uteži po enačbi (5.20) h = 0.3 Uteži po enačbi (5.21) h = 0.1 Uteži po enačbi (5.19) h = 0.3 Uteži po enačbi (5.20) h = 0.3 Uteži po enačbi (5.21) h = 0.1 SLIKA 5.3: Prostorsko utežena regresija z različnimi funkcijami w(d) S ciljem, da bi ocenili, kateri način prostorsko utežene regresije je najboljši, smo opravili postopek navzkrižne validacije. Za vseh 359 točk, v katerih poznamo vrednosti slučajnega polja (koncentracije kroma) smo izračunali prostorsko uteženo regresijo, pri čemer pa nismo upoštevali vrednosti v točki, za katero smo opravili regresijo. Za vsako od 359 točk smo torej uporabili podatke na preostalih 358 točkah. Tako pridobljene ocene vrednosti slučajnega polja smo primerjali z izmerjenimi in po enačbi (5.22) izračunali vrednost navzkrižne validacije CV . Rezultate za vse tri predlagane enačbe ((5.19)– (5.21)) in različne vrednosti h prikazujemo na sliki 5.4. Vidimo lahko, da je dala enačba (5.21) najboljše rezultate, saj so tu razlike med izmerjenimi in napovedanimi vrednostmi slučajnega polja najmanjše. S slike tudi opazimo, da je za enačbi (5.19) in (5.20) optimalna razdalja h približno 0.3 km, medtem ko je za enačbo (5.21) optimalna razdalja h približno 0.1 km. Ti vrednosti smo uporabili tudi za prikaz prostorske regresije na sliki 5.3. 68 5 Prostorsko utežena regresija 45 000 CV 40 000 35 000 enačba (5.19) 30 000 enačba (5.20) enačba (5.21) 25 0000 1 2 3 4 h SLIKA 5.4: Določitev optimalne vrednosti h z uporabo navzkrižne validacije 6 Prostorski vzorci Pri analizi prostorskih podatkov so lege točk, kjer merimo vrednosti slučajne funkcije oziroma slučajnega polja, zelo pomemben podatek. Razpored točk je lahko rezultat načrta poskusa ali pa je slučajen. Zaradi različnih dejavnikov so lahko točke enakomerno razporejene po prostoru ali pa so v gručah oziroma klastrih. Pri razporedu točk po prostoru lahko govorimo o treh značilnih primerih: • slučajno enakomerno porazdeljene točke v prostoru, • slučajno porazdeljene točke v gručah (v nadaljevanju v gručah), • deterministično enakomerno razporejene točke v prostoru. (a) (b) (c) SLIKA 6.1: Trije tipi razporeditve točk: (a) slučajno enakomerno, (b) v gručah in (c) enakomerno Osnovna naloga tega področja prostorske statistike je, da znamo ločiti med temi tremi primeri oziroma da lahko s statističnimi preizkusi in izbranim tveganjem α zavržemo ali ne zavržemo določene domneve o porazdeljenosti točk v slučajnem polju in zapišemo odločitve. Pri preizkušanju domnev bomo za ničelno domnevo vedno izbrali izjavo, da so točke porazdeljene slučajno enakomerno. Oblika alternativne domneve pa je odvisna od tega, ali opravljamo enostranski ali dvostranski preizkus. Poznamo torej tri možnosti ničelne in alternativne domneve: H0: Točke so razporejene slučajno enakomerno. H1: Točke niso razporejene slučajno enakomerno. ali 70 6 Prostorski vzorci H0: Točke so razporejene slučajno enakomerno. H1: Točke so razporejene v gručah. ali H0: Točke so razporejene slučajno enakomerno. H1: Točke so razporejene deterministično enakomerno. Analizo prostorskih vzorcev lahko naredimo z metodo kvadratov ali pa z metodo razdalj oziroma metodo najbližjega soseda. 6.1 Metoda kvadratov Pri metodi kvadratov na območje analize postavimo kvadrate ali druge like ter preštejemo, koliko točk je znotraj posameznih likov. Na sliki 6.2 prikazujemo različne možnosti postavitev likov različnih oblik. Lahko se odločimo za povsem slučajno postavitev likov ali pa like postavimo tako, da z njimi prekrijemo celotno območje. Ne glede na obliko lika, bomo v okviru razlag te metode vse like imenovali ,,kvadrati“. SLIKA 6.2: Različne oblike likov, ki jih lahko uporabimo pri metodi kvadratov Število točk v posameznem kvadratu označimo z Yi, i = 1, . . . , n, kjer je n število kvadratov, ki jih 6.1 Metoda kvadratov 71 uporabimo v analizi. Ob predpostavki, da so točke razporejene slučajno enakomerno, je verjetnost, da je točka znotraj kvadrata, enaka razmerju ploščin kvadrata Ak in celotnega območja A: Ak p = , (6.1) A število točk znotraj kvadrata pa je slučajna spremenljivka, porazdeljena po Poissonovi porazdelitvi. Verjetnostna funkcija Poissonove porazdelitve je: νy e−ν pY (y) = , y = 0, 1, 2, . . . , (6.2) y! kjer parameter ν predstavlja pričakovano število točk v posameznem kvadratu in je enako ν = p N . Z N smo označili število vseh točk na obravnavanem območju. Zanimiva lastnost Poissonove porazdelitve je, da sta srednja vrednost in varianca enaki parametru porazdelitve ν: E[Y ] = var[Y ] = ν. (6.3) Na osnovi te lastnosti lahko pripravimo eno od možnih statistik za preizkus domneve o prostorski ure-jenosti: S∗2 H∗ = Y ∼ χ∗2 ¯ Y n−1, (6.4) kjer smo s S∗2 označili varianco vzorca, z ¯ Y pa njegovo povprečje. Statistika H∗ se porazdeljuje po Y normirani porazdelitvi χ∗2, ki se od običajne porazdelitve χ2 loči po tem, da ima srednjo vrednost enako ena, oziroma da je H∗ le (n − 1)-krat manjša od H: H H∗ = −→ F n − 1 χ∗2 (h) = Fχ2 ((n − 1)h). (6.5) V zadnji enačbi smo upoštevali enačbo za določitev porazdelitve funkcije slučajne spremenljivke za primer monotono naraščajoče funkcije, kar zveza v (6.5) za n > 1 seveda je. Enačbo za določitev porazdelitvene funkcije slučajne spremenljivke Y , ki je monotono naraščajoča funkcija X, namreč zapišemo takole: Y = g(X) −→ FY (y) = FX (g−1(y)). (6.6) Poleg statistike H∗, ki jo določimo z enačbo (6.4), so različni avtorji uporabili tudi druge statistike. Te so na primer: S∗2 ICS = Y − 1 ¯ (David, Moore, 1954), Y ˙ S∗2 Y = ¯ Y + Y − 1 ¯ (Lloyd,1967), Y ˙ Y (6.7) IP = ¯ (Lloyd,1967), Y n P Yi(Yi − 1) I i=1 δ = ¯ (Morisita, 1959). Y (n ¯ Y − 1) 72 6 Prostorski vzorci Porazdelitve vseh teh statistik niso splošno znane. Eden od načinov za uporabo le-teh je uporaba simulacij. Porazdelitev statistike ICS je sicer podobna porazdelitvi χ∗2, le da je zamaknjena za ena proti levi. Za porazdelitev statistike ˙ Y lahko za nekoliko večje število kvadratov, na primer, če je število kvadratov večje od 80, uporabimo normalno porazdelitev s srednjo vrednostjo ν in standardnim odklonom pν/n + 2/(n − 1). Pri metodi kvadratov lahko uporabimo eno od statistik v enačbah (6.4) ali (6.7), lahko pa uporabimo preizkus domneve o porazdelitvi. Če vemo, da se število točk v posameznem kvadratu ob predpostavki, da so točke razporejene slučajno enakomerno, porazdeljuje po Poissonovi porazdelitvi, lahko ničelno domnevo preizkusimo s testom Poissonove porazdelitve. To lahko naredimo s preizkusom χ2, s preizkusom Kolmogorova in Smirnova ali kakšno drugo metodo za preizkus domnev o porazdelitvi. Pri preizkusu χ2 preštejemo, koliko kvadratov vsebuje določeno število točk. Glede na pričakovano število točk znotraj kvadrata in število kvadratov se odločimo, na koliko razredov razdelimo podatke in kaj predstavljajo ti razredi. Včasih je primerno, da prvi razred predstavlja kvadrate, ki nimajo nobene točke, drugi kvadrate z eno točko in tako naprej. Zadnji K-ti razred pa predstavlja vse tiste kvadrate, ki vsebujejo K − 1 ali več točk. Če je pričakovano število točk v kvadratu večje, je smiselno, da prvi razred predstavlja kvadrate, ki vključujejo manj kot določeno število točk. V preglednici 6.1 prikazujemo dva značilna primera razporeditev razredov. PREGLEDNICA 6.1: Dva primera razporeditve razredov glede na število točk v kvadratu Razred Število točk Razred Število točk 1 0 1 0–3 2 1 2 4 3 2 3 5 4 3 4 6 5 4 5 7 6 5 6 8–9 7 6 ali več 7 10 ali več Za vsak razred izračunamo teoretična števila kvadratov, ki vsebujejo določeno število točk (teoretične velikosti razredov) ni: ni = P [število točk v kvadratu ustreza razredu i] n, (6.8) kjer je n število vseh kvadratov. Empirične velikosti razredov ˆ ni določimo tako, da preštejemo kvadrate, ki ustrezajo razredu i. Nato lahko izračunamo statistiko H po naslednji enačbi: K X (ni − ˆ ni)2 H = ∼ χ2 n K−1, (6.9) i i=1 6.2 Primeri – metoda kvadratov 73 kjer smo označili, da se statistika H porazdeljuje po porazdelitvi χ2 s K − 1 prostostnimi stopnjami. Pri tem preizkusu ne moremo ločiti med razporeditvijo v gruče in urejeno deterministično razporeditvijo. Ugotovimo lahko le, ali razporeditev točk statistično značilno odstopa od slučajne enakomerne razporeditve. Kritično območje je H > χ2 . Če je vrednost statistike v kritičnem območju, K−1,1−α moramo ničelno domnevo zavrniti. 6.2 Primeri – metoda kvadratov Analizirajmo tri različne primere razporeditve točk: slučajno enakomerno (primer A), v gručah (primer B) in deterministično enakomerno (primer C). Postavimo ničelno in alternativno domnevo ter izberimo tveganje α, pri katerem bi še zavrnili ničelno domnevo: H0: Točke so razporejene slučajno enakomerno. H1: Točke niso razporejene slučajno enakomerno. α = 0.05 Razporeditev točk prikazujemo na slikah 6.3 in 6.4. Kvadratno območje razdelimo na 16 kvadratov in preštejemo število točk v posameznih kvadratih. Skupno število točk N je enako 100 za vse tri primere. Na slikah so označeni tudi kvadrati s številkami od 1 do 16. 4 4. 8. 12. 16. Primer A: Slučajna razporeditev 3 3. 7. 11. 15. y 2 2. 6. 10. 14. 1 1. 5. 9. 13. 00 1 2 3 4 x SLIKA 6.3: Enakomerna slučajna razporeditev točk – primer A 74 6 Prostorski vzorci 4 4 4. 8. 12. 16. 4. 8. 12. 16. Primer B: Razporeditev v gručah Primer C: Enakomerna razporeditev 3 3 3. 7. 11. 15. 3. 7. 11. 15. y y 2 2 2. 6. 10. 14. 2. 6. 10. 14. 1 1 1. 5. 9. 13. 1. 5. 9. 13. 0 0 0 1 2 3 x 4 0 1 2 3 4 x SLIKA 6.4: Različne razporeditve točk – primera B in C V preglednici 6.2 je zapisano, koliko točk je v posameznem kvadratu. Iz teh podatkov lahko izračunamo povprečje in varianco ter določimo statistiko H∗ po enačbi (6.4). PREGLEDNICA 6.2: Število točk v posameznih kvadratih Kvadrat 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Primer A 7 7 3 4 4 12 4 5 4 7 7 9 5 9 6 7 Primer B 3 3 4 2 7 10 9 3 3 10 13 1 5 5 10 0 Primer C 8 7 7 8 6 4 5 5 7 8 7 8 4 6 4 6 V preglednici 6.3 prikazujemo vrednosti statistike z mejami kritičnega območja. V zadnjem stolpcu so vrednosti dejanskega tveganja ob zavrnitvi ničelne domneve. V primeru, da je statistika manjša od spodnje meje H∗ < H1 ali večja od zgornje meje H∗ > H2, ničelno domnevo zavrnemo. V teh primerih je dejansko tveganje αdej < α = 0.05. PREGLEDNICA 6.3: Izračun statistik in kritičnih območij ¯ Y S∗2 H∗ H Y 1 H2 αdej Primer A 6.25 6.42 1.01 0.417 1.83 0.874 Primer B 5.50 14.8 2.69 0.417 1.83 0.000799 Primer C 6.25 0.77 0.35 0.417 1.83 0.021 6.2 Primeri – metoda kvadratov 75 V naših primerih lahko ugotovimo naslednje: v primeru A statistika ni v kritičnem območju, dejansko tveganje je večje od 5 %, saj znaša kar 87.4 %, ničelne domneve ne moremo zavrniti in s tveganjem α ne moremo trditi, da točke niso porazdeljene slučajno enakomerno. V primerih B in C je enkrat statistika večja od zgornje meje, drugič pa manjša od spodnje. V obeh primerih je dejansko tveganje manjše od 5 %, ničelno domnevo moramo zavrniti in s tveganjem α, manjšim od 5 %, lahko trdimo, da točke niso porazdeljene slučajno enakomerno. Analizirajmo še porazdelitev točk, v katerih so bile opravljene meritve v območju Jura. Naredili bomo enostranski preizkus, zato je alternativna domneva drugačna, kot v prejšnjem primeru: H0: Točke so razporejene slučajno enakomerno. H1: Točke so zbrane v gručah. α = 0.05 Nalogo ponovimo z dvema različnima postavitvama kvadratov. V prvem primeru uporabimo kvadrate velikosti 1 km × 1 km tako, da je bilo pokrito območje točk, v drugem primeru pa uporabimo manjše kvadrate velikosti 0.5 km × 0.5 km, razporejeni pa so tako, kot kaže slika 6.5 (desno). 6 6 5 5 y [km] y [km] 4 4 3 3 2 2 1 1 0 0 0 1 2 3 4 5 0 1 2 3 4 5 x [km] x [km] SLIKA 6.5: Dva načina postavitve kvadratov V preglednici 6.4 prikazujemo izračun statistike, mej kritičnega območja in dejanskega tveganja. Ugotovimo, da postavitev kvadratov odločilno vpliva na končen rezultat in odločitev o tipu razporeditve točk. Vidimo, da je ob prvi postavitvi kvadratov tveganje za zavrnitev ničelne domneve zelo velika. Zato ne moremo trditi, da razporeditev točk ni slučajna enakomerna. Šele ob uporabi druge postavitve manjših 76 6 Prostorski vzorci kvadratov metoda prepozna, da so točke združene v gruče. V tem primeru je statistika v kritičnem območju, saj je H∗ = 2.10 > H2 = 1.48. Dejansko tveganje je seveda manjše od dovoljenega αdej = 0.0006 ≪ α = 0.05. Zato lahko s tveganjem, ki je precej manjše od 5 %, trdimo, da se točke združujejo v gruče. PREGLEDNICA 6.4: Izračun statistik in kritičnih območij ¯ Y S∗2 H∗ H Y 2 αdej Postavitev A 24.7 25.9 1.05 1.72 0.400 Postavitev B 7.21 15.1 2.10 1.48 0.0006 Pri razporedu A število kvadratov ne zadošča za uporabo preizkusa χ2, zato ga bomo opravili le za razpored B. Odločimo se, da bodo v prvem razredu kvadrati, ki vsebujejo od nič do štiri točke, v drugem kvadrati s pet točkami, v tretjem kvadrati s šest točkami in tako naprej do šestega razreda s kvadrati, ki vsebujejo devet točk. V sedmem, zadnjem razredu so kvadrati, ki imajo več kot devet točk. V preglednici 6.5 prikazujemo razporeditev v razrede, teoretične velikosti razredov ni in empirične velikosti razredov ˆ ni. PREGLEDNICA 6.5: Izračun statistik in kritičnih območij razred i razporeditev ni ˆ ni 1 0–4 4.493 16 2 5 3.483 4 3 6 4.184 4 4 7 4.308 0 5 8 3.881 2 6 9 3.108 0 7 10 ali več 5.542 3 Postavimo ničelno in alternativno domnevo ter izberimo tveganje α, pri katerem bi še zavrnili ničelno domnevo: H0: Točke so razporejene slučajno enakomerno. H1: Točke niso razporejene slučajno enakomerno. α = 0.05 Statistiko H izračunamo po enačbi (6.9) in dobimo: (4.493 − 16)2 (3.483 − 4)2 (5.542 − 3)2 H = + + · · · + = 39.05. (6.10) 4.493 3.483 5.542 6.3 Pomanjkljivosti metode kvadratov 77 Meja kritičnega območja je χ2 = 12.59 6,0.95 . Ker je vrednost statistike bistveno višja od meje kritičnega območja H = 39.05 > 12.59, dejansko tveganje pa je zanemarljivo αdej = 7.01 · 10−7, moramo ničelno domnevo zavrniti in s tveganjem, manjšim od 5 %, trdimo, da točke niso razporejene slučajno enakomerno. 6.3 Pomanjkljivosti metode kvadratov Zadnji primer, ko smo zaradi različne razporeditve in velikosti kvadratov prišli do povsem nasprotujočih si rezultatov, kaže na pomanjkljivosti metode kvadratov. Te lahko opredelimo v naslednjih točkah: • Rezultat je lahko odvisen od velikosti in razporeditve kvadratov. Zato je pri uporabi metode kvadratov primerno, da se preizkus ponovi pri različnem številu, velikosti in razporeditvi kvadratov. • Ker statistiko izračunamo iz števila točk v kvadratih, dejansko primerjamo le povprečne gostote točk v kvadratih, ne pa razporeditev ene točke glede na druge. • Pri analizi izgubimo informacijo o razporeditvi točk znotraj kvadrata. Zato gruč, ki so manjše od kvadratov, ta metoda ne zazna. Da bi premostili pomanjkljivosti metode kvadratov, lahko uporabimo alternativno metodo najbližjega soseda. 6.4 Metoda najbližjega soseda Pri metodi najbližjega soseda za vsako točko izračunamo razdaljo d do najbližje sosednje točke: n di = min dij, i ̸= j. (6.11) j=1 Na sliki 6.6 so prikazane točke z enakomerno razporeditvijo in razporeditvijo v gruče. Prikazane so tudi razdalje od posamezne točke do najbližje sosede. Že površna primerjava odkrije, da so razdalje do najbližjih sosedov večinoma krajše v primeru, ko so točke razporejene v gručah. 78 6 Prostorski vzorci 1.0 1.0 0.8 0.8 y y 0.6 0.6 0.4 0.4 0.2 0.2 Enakomerna slučajna razporeditev Razporeditev v gručah 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 x x SLIKA 6.6: Najbližji sosedi za dve različni razporeditvi točk V primeru, da je razporeditev točk slučajna enakomerna, porazdelitev razdalje do najbližjega soseda lahko zapišemo z naslednjo gostoto verjetnosti: fD(d) = 2π λ d e−π λ d2, za d > 0, (6.12) kjer n λ = (6.13) A predstavlja gostoto točk (n je število točk, A pa ploščina analiziranega območja). Graf te gostote verjetnosti prikazujemo na sliki 6.7. Zapišimo še srednjo vrednost in varianco slučajne spremenljivke D: r r 1 A 1 1 mD = E[D] = = , 2 n 2 λ (6.14) A 1 1 1 1 1 σ2 − − D = var[D] = = . n2 π 4 nλ π 4 6.4 Metoda najbližjega soseda 79 1.5 fD(d) 1.0 0.5 0.00.0 0.5 1.0 1.5 2.0 d SLIKA 6.7: Gostota verjetnosti razdalje do najbližjega soseda za λ = 1 Pri metodi najbližjega soseda lahko izberemo tudi alternativno pot, da namesto iskanja najbližjega soseda vsake točke določamo razdaljo od navidezne točke do njej najbližje dejanske točke. Pri tem predpostavimo, da je navidezna točka kjerkoli na območju, ki ga analiziramo. To pomeni, da je porazdelitev lege navidezne točke slučajna enakomerna. Na sliki 6.8 prikazujemo najkrajše razdalje od navideznih točk (obarvanih sivo) do dejanskih točk za dve različni razporeditvi točk. Tudi tu lahko že na prvi pogled opazimo razlike med razdaljami do točk, če so dejanske točke razporejene enakomerno ali v gručah. 1.0 1.0 0.8 0.8 y y 0.6 0.6 0.4 0.4 0.2 0.2 Enakomerna slučajna razporeditev Razporeditev v gručah 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 x SLIKA 6.8: Najkrajša razdalja od navideznih do analiziranih točk za dve razporeditvi dejanskih točk Porazdelitev razdalje X od navidezne do najbližje dejanske točke je enaka porazdelitvi razdalje D do najbližjega soseda za primer slučajne enakomerne razporeditve, saj navidezna točka, od katere merimo razdaljo, predstavlja le še eno točko med slučajno enakomerno razporejenimi točkami. Edina, skoraj zanemarljiva razlika je v izračunu gostote točk, ki je v tem primeru enaka n + 1 λ = . (6.15) A 80 6 Prostorski vzorci Možnosti za določitev ustrezne statistike, na osnovi katere bomo lahko izpeljali statistični preizkus, je več. Temeljijo tako na slučajni spremenljivki D, to je razdalji med najbližjimi sosedami, kot tudi na slučajni spremenljivki X, torej razdalji med navidezno in njej najbližjo dejansko točko. Zapišimo nekaj različnih statistik z ustreznimi porazdelitvami: √ n 2 λ X 4 − π U = di ∼ N 1, (Clark, Evens, 1954), (6.16) n nπ i=1 n X H = 2πλ d2 ∼ i χ22n (Skellam, 1952), (6.17) i=1 n πλ X 1 U = x2 ∼ N 1, (Pielou, 1959), (6.18) n i n i=1 n n P x2i porazdelitve ne poznamo, uporabiti S = i=1 ∼ (6.19) n 2 moramo simulacije (Eberhardt, 1967), P xi i=1 n ! n ! 12n 1 X X H = n ln x2 − ln x2 ∼ χ2 7n + 1 n i i n−1 (Pollard, 1971), (6.20) i=1 i=1 n 1 X x2 1 1 U = i ∼ N , (Byth, Ripley, 1980), (6.21) n x2 + d2 2 12n i=1 i i n P x2i F = i=1 ∼ F n 2n,2n (Hopkins, 1954). (6.22) P d2i i=1 Nekateri avtorji pa so pri definiciji statistike uporabili še slučajno spremenljivko X2, ki predstavlja razdaljo med navidezno in drugo najbližjo dejansko točko. n 1 X x2 1 1 U = i ∼ N , (Holgate, 1965), (6.23) n x2 2 12n i=1 2i n P x2i B = i=1 ∼ β n n,n (Holgate, 1965), (6.24) P x22i i=1 kjer xi predstavlja razdaljo med navidezno točko in njej najbližjo točko, x2i pa razdaljo med navidezno točko in njej drugo najbližjo točko. 6.5 Primeri – metoda najbližjega soseda 81 6.5 Primeri – metoda najbližjega soseda Analizirajmo enake primere kot pri metodi kvadratov v podpoglavju Primeri - metoda kvadratov. Obravnavamo tri različne primere razporeditve točk na kvadratnem območju velikosti 4 × 4: slučajno enakomerno (primer A), v gručah (primer B) in deterministično enakomerno (primer C). Postavimo ničelno in alternativno domnevo ter izberimo tveganje α, pri katerem bi še zavrnili ničelno domnevo: H0: Točke so razporejene slučajno enakomerno. H1: Točke niso razporejene slučajno enakomerno. α = 0.05 Pripravili smo vzorce s stotimi točkami za različne razporeditve dejanskih točk. Nato smo pripravili še vzorec stotih navideznih točk, ki so razporejene slučajno enakomerno. Iz teh vzorcev smo lahko določili razdalje di, xi in x2i, iz katerih lahko izračunamo statistike po enačbah (6.16)–(6.24). 82 6 Prostorski vzorci Različne razporeditve stotih točk grafično prikazujemo na sliki 6.9, kjer so označene najkrajše razdalje med dejanskimi točkami – vzorec di. 4 4 Primer A: Enakomerna slučajna 3 3 y y 2 2 1 1 Primer B: Razporeditev v gručah 0 0 0 1 2 3 4 x 0 1 2 3 4 x 4 Primer C: Urejena razporeditev 3 y 2 1 00 1 2 3 4 x SLIKA 6.9: Najkrajše razdalje med točkami za tri razporeditve dejanskih točk 6.5 Primeri – metoda najbližjega soseda 83 Na sliki 6.10 prikazujemo podobno razporejene točke kot na sliki 6.9, poleg njih pa še enako veliko skupino slučajno enakomerno razporejenih navideznih točk, ki so označene s svetlejšo sivo barvo. Tu označujemo še razdalje od navideznih točk do najbližje dejanske točke z zelenimi črtami – vzorec xi in razdalje od navideznih do druge najbližje dejanske točke z rdečimi črtami – vzorec x2i. 4 4 y y 3 3 2 2 1 1 Primer A: Enakomerna slučajna Primer B: Razporeditev v gručah 00 1 2 3 4 0 x 0 1 2 3 4 x 4 y 3 2 1 Primer C: Urejena razporeditev 00 1 2 3 4 x SLIKA 6.10: Najkrajše razdalje med navideznimi in dejanskimi točkami za tri razporeditve točk V preglednici 6.6 prikazujemo izračun statistik, mej kritičnega območja in dejanskega tveganja za tri različne tipe razporeditve točk in za osem različnih načinov določitve statistike. Vidimo lahko, da je v primeru A, ko so dejanske točke porazdeljene slučajno enakomerno, statistika vedno v srednjem območju, ki nam ne narekuje zavrnitve ničelne domneve. Zato v tem primeru ničelne domneve ne zavrnemo in zaključimo, da s tveganjem α = 5 % ne moremo trditi, da točke niso razporejene slučajno ena- 84 6 Prostorski vzorci komerno. V tem primeru je za vse različne statistike dejansko tveganje večje od izbranega. V primeru B vidimo, da je statistika v kritičnem območju vedno, razen pri uporabi statistike (6.21). V primerih, ko je statistika v kritičnem območju in je s tem tudi tveganje manjše od izbranega tveganja 5 %, ničelno domnevo zavrnemo in zaključimo, da s tveganjem, manjšim od 5 %, lahko trdimo, da točke niso razporejene slučajno enakomerno. V primeru C, ko so točke razporejene deterministično enakomerno, vsi preizkusi kažejo na to, da moramo ničelno domnevo zavrniti, saj je statistika v kritičnem območju, dejansko tveganje pa nižje od izbranega. Tudi v tem primeru pridemo do enakega zaključka, da s tveganjem, manjšim od 5 %, lahko trdimo, da točke niso razporejene slučajno enakomerno. PREGLEDNICA 6.6: Izračun statistik, kritičnih območij in dejanskega tveganja za tri razporeditve točk Enačba spodnja zgornja Primer A Primer B Primer C meja meja statistika tveganje statistika tveganje statistika tveganje (6.16) 0.898 1.102 1.022 0.668 0.789 5.61 · 10−5 1.827 0. (6.17) 162.7 241.1 222.3 0.268 160.2 0.0351 535. 0. (6.18) 0.804 1.196 0.959 0.688 4.093 0. 0.639 0.000303 (6.20) 73.36 128.4 88.98 0.490 238.5 3.23 · 10−13 63.77 0.00460 (6.21) 0.443 0.557 0.493 0.818 0.541 0.156 0.404 0.000904 (6.22) 0.431 0.569 0.456 0.215 0.684 7.78 · 10−8 0.331 8.15 · 10−7 (6.23) 0.443 0.557 0.499 0.976 0.651 1.79 · 10−7 0.190 0. (6.24) 0.757 1.320 0.863 0.301 5.111 0. 0.239 0. Naredimo še preizkus domneve o razporeditvi točk, v kateri so bile opravljene meritve na območju Jura. Naredimo enak preizkus kot za prejšnji primer, ko smo obravnavali porazdelitev točk v kvadratu 4 × 4. Število točk je enako 359, prav toliko točk tudi generiramo, in sicer tako, da so razporejene slučajno enakomerno. Gostota točk je λ = 23.3 točk/km2. Razporeditev točk, kjer so bile opravljene meritve, in enakomerno porazdeljenih navideznih točk ter razdalje do najbližjih točk so prikazane na sliki 6.11. Na tej sliki so s sivo označene razdalje med najbližjimi sosednjimi točkami, kjer so bile opravljene meritve, z zeleno pa razdalje od navideznih do najbližjih točk meritev. H0: Točke so razporejene slučajno enakomerno. H1: Točke niso razporejene slučajno enakomerno. α = 0.05 6.5 Primeri – metoda najbližjega soseda 85 6 Točke meritev 5 Navidezne točke y [km] 4 3 2 1 00 1 2 3 4 5 x [km] SLIKA 6.11: Razporeditev točk meritev in navideznih točk za primer Jura V preglednici 6.7 prikazujemo izračune za preizkušanje domneve o razporeditvi merskih točk za primer Jura. Ugotovimo lahko, da metode različnih avtorjev pripeljejo do različnih zaključkov. V primeru uporabe enačb (6.16), (6.17) ter (6.20)–(6.22) je tveganje manjše od izbrane vrednosti α = 5 %, zato lahko ničelno domnevo zavrnemo in izjavimo, da s tveganjem, manjšim od 5 %, lahko trdimo, da točke niso razporejene slučajno enakomerno. Če pa računamo po enačbi (6.18), ničelne domneve ne moremo zavrniti. V tem primeru izjavimo, da ob tveganju 5 % ne moremo trditi, da točke niso razporejene slučajno enakomerno. Poudariti moramo, da so statistike, pri katerih uporabimo navidezne točke, odvisne od slučajno generiranega vzorca navideznih točk. Zato ponovitve teh izračunov lahko vodijo do različnih rezultatov in zato tudi do različnih zaključkov. 86 6 Prostorski vzorci PREGLEDNICA 6.7: Preizkušanje domneve o razporeditvi točk v primeru Jura Enačba spodnja meja zgornja meja statistika tveganje (6.16) 0.946 1.054 1.209 3.62 · 10−14 (6.17) 645.6 794.1 1358. 0. (6.18) 0.897 1.103 1.084 0.114 (6.20) 307.5 412.3 302.6 0.0305 (6.21) 0.470 0.530 0.550 0.00115 (6.22) 0.864 1.158 0.551 2.32 · 10−15 7 Metoda Monte Carlo Integriranje je zelo pomembna in uporabna operacija na mnogih področjih tehnike, tudi v prostorski statistiki. Ker sta v mnogih primerih tako integrant kot integracijsko območje taka, da le za posamezne preproste primere ta integral lahko izračunamo analitično, moramo uporabiti numerične metode. Metoda Monte Carlo je ena od metod numeričnega integriranja, ki temelji na generiranju vzorca slučajnih spremenljivk in vektorjev. Osnovni postopek metode Monte Carlo je preprost: 1. Definiramo problem, definiramo osnovne lastnosti slučajnih spremenljivk, ki jih uporabimo pri reševanju problema. Določimo število simulacij (število ponovitev). 2. Začetek zanke simulacij (ponavljanje poskusov). 2.1 Generiramo vzorec slučajnih spremenljivk. 2.2 Opravimo analizo, s katero rešimo osnovni problem v okviru ene simulacije. 3. Konec zanke simulacij. 4. Opravimo končno analizo ter numerični in grafični prikaz rezultatov. Vzorec slučajne spremenljivke tu razumemo kot vrednosti, ki jih slučajna spremenljivka zavzame ob ponavljanju poskusa. Vzemimo, da želimo s to metodo določiti verjetnost, da je rezultat meta poštene kocke manjši ali enak dve. Rezultat seveda poznamo in je enak 1/3: 1. Pri definiciji problema zapišemo, kaj želimo določiti in da je slučajna spremenljivka diskretna in enakomerno porazdeljena od 1 do vključno 6. Določimo število simulacij, na primer 10000. Postavimo števec ugodnih izidov na nič. 2. Začetek zanke simulacij (ponavljanje poskusov). 2.1 Generiramo številko od ena do šest po diskretni enakomerni porazdelitvi. 2.2 Ugotovimo, ali je generirana številka manjša ali enaka dve. Če je, povečamo števec ugodnih izidov za ena. 3. Konec zanke simulacij. 4. Oceno iskane verjetnosti dobimo tako, da število ugodnih izidov delimo s številom vseh simulacij. 88 7 Metoda Monte Carlo Ime metode Monte Carlo se je prvič pojavilo med drugo svetovno vojno, ko so znanstveniki von Neu-man, Ulam, Fermi in drugi v laboratorijih v Los Alamosu tudi s pomočjo simulacij skonstruirali prvo atomsko bombo. Sicer pa sta ideja o generiranju slučajnih spremenljivk in njihova uporaba pri reševanju problemov bistveno starejši. Prvi dobro dokumentiran primer je opravil Comte de Buffon leta 1777, ko je v problemu, poznanem kot Buffonova igla, s simulacijami poskušal določiti vrednost števila π. Pozneje so generiranje slučajnih števil uporabili tudi Lord Kelvin (1901), Gosset (Student) (1908), Fermi (1930) in drugi (Kalos et al., 1986, Rubinstein, 1981). Z razvojem boljših, hitrejših računalnikov postaja metoda Monte Carlo uporabna za vse širši krog ljudi. Trenutno, v jeseni 2023, lahko v sistemu Web of Science najdemo več kot 143000 člankov, ki vključujejo področje metode Monte Carlo. Pri primeru, ki ga poznamo kot Buffonova igla, na površino, na kateri so narisane vzporedne črte z medsebojno razdaljo d, postopoma mečemo igle. Comte de Buffon je pokazal, da je verjetnost, da igla dolžine l seka eno izmed vzporednih črt, enaka: 2 l 2 l P [igla seka črto] = → π = . π d P [igla seka črto] d Če torej s poskušanjem ocenimo verjetnost, da igla seka črto, lahko po zgornji enačbi ocenimo število π. Animacija takega postopka določitve števila π je prikazana na sliki 7.1. n = 300 π = 3.226 SLIKA 7.1: Animacija določitve števila π 7.1 Generiranje vzorca slučajne spremenljivke 89 7.1 Generiranje vzorca slučajne spremenljivke Generiranje vzorca slučajne spremenljivke X temelji na vzorcu enakomerno porazdeljene slučajne spremenljivke U s porazdelitveno funkcijo FU (u) = u (za 0 ≤ u ≤ 1). Vzorca povsem slučajne spremenljivke ni lahko generirati. Zato z računalniki običajno generiramo vzorce tako imenovanih psevdoslučajnih števil, ki so pravzaprav zaporedja determinističnih števil z zelo veliko periodo ponovitve. Ta zaporedja imajo enake lastnosti kot zaporedja povsem slučajnih števil. Skoraj vsi računalniški programi, namenjeni računanju (na primer prevajalniki: FORTRAN, C++, PASCAL, BASIC in drugi programi: MATLAB, MATHEMATICA, EXCEL, R, PYTHON), vključujejo tudi ,,generatorje slučajnih števil“, s katerimi generiramo zaporedja psevdoslučajnih števil, ki ustrezajo vzorcu slučajne spremenljivke U , ki je porazdeljena enakomerno od 0 do 1. Vzorec poljubno porazdeljene slučajne spremenljivke ali vektorja lahko generiramo z različnimi metodami: z inverzno metodo, metodo sprejema in zavrnitve, polarno metodo, metodo trakov, z mrežno metodo in drugimi. Najpogosteje uporabljamo inverzno metodo. 7.1.1 Inverzna metoda generiranja vzorca Vzemimo, da imamo enako velika vzorca dveh slučajnih spremenljivk. Prva slučajna spremenljivka, X, je porazdeljena po porazdelitvi FX (x), druga, U , pa je enakomerno porazdeljena od 0 do 1: FU (u) = u (za 0 ≤ u ≤ 1). Elemente v obeh vzorcih razvrstimo po velikosti in predpostavimo, da sta verjetnosti, da sta slučajni spremenljivki X in U manjši od i-tega elementa ustreznega vzorca, enaki, ne glede na porazdelitev. To predpostavko v enačbi za vsak i zapišemo takole: P [X ≤ xi] = FX (xi) = P [U ≤ ui] = FU (ui) = ui. Iz zadnje enačbe lahko izpeljemo izraz za določitev elementa vzorca slučajne spremenljivke: FX (xi) = ui −→ xi = F −1(u X i). (7.1) Če torej poznamo porazdelitveno funkcijo slučajne spremenljivke oziroma njeno inverzno funkcijo in imamo vzorec (ui, i = 1, . . . , n) enakomerno porazdeljene slučajne spremenljivke U , lahko po zadnji enačbi določimo vzorec slučajne spremenljivke X, kjer je n velikost vzorca. Inverzno metodo na-zorno prikažemo tudi grafično na sliki 7.2. Z modro je označena enakomerno porazdeljena slučajna spremenljivka U s pripadajočim vzorcem. Vsakemu elementu vzorca slučajne spremenljivke U ustreza element vzorca slučajne spremenljivke X. Na zgornjem delu slike metodo opišemo na primeru treh generiranih vrednosti, na spodnjem delu pa je prikazan nekoliko večji vzorec. Na spodnji sliki lahko opazimo, da so točke v okolici srednje vrednosti mX = 6 bolj zgoščene kot na obeh repih, kjer je točk manj. 90 7 Metoda Monte Carlo 1.2 Fu(u) Fx(x) 1.0 Enakomerna 0.8 porazdelitev 0.6 0.4 Normalna porazdelitev 0.2 0.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 12 14 u x 1.2 Fu(u) Fx(x) 1.0 Enakomerna porazdelitev 0.8 0.6 0.4 Normalna 0.2 porazdelitev 0.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 12 14 u x SLIKA 7.2: Inverzna metoda generiranja vzorca V preglednici 7.1 podajamo nekaj značilnih primerov uporabe inverzne metode za generiranje vzorcev slučajnih spremenljivk, katerih porazdelitvene funkcije poznamo v zaključeni obliki. PREGLEDNICA 7.1: Uporaba inverzne metode Porazdelitev Način generiranja Dogodek, ki se zgodi z verjetnostjo p Dogodek se zgodi, če je ui manjši od p Enakomerna zvezna X od a do b xi = a + ui (b − a) Eksponentna porazdelitev T ti = − 1 ln(1 − u λ i) Gumbelova porazdelitev maksimuma Y yi = u − 1 ln(− ln u α i) Gumbelova porazdelitev minimuma Z zi = u + 1 ln(− ln(1 − u α i)) Fréchetova porazdelitev maksimuma Y yi = u (− ln ui)− 1k Fréchetova porazdelitev minimuma Z zi = u (− ln(1 − ui))− 1k 1 Weibullova porazdelitev maksimuma Y yi = ω − (ω − u) (− ln ui) k 1 Weibullova porazdelitev minimuma Z zi = ε + (u − ε) (ln(1 − ui)) k 7.1 Generiranje vzorca slučajne spremenljivke 91 7.1.2 Generiranje vzorca normalno porazdeljenih slučajnih spremenljivk Porazdelitvene funkcije normalne porazdelitve v zaključeni obliki ne poznamo, prav tako ne poznamo njene inverzne funkcije. Zato za generiranje vzorca normalno porazdeljene slučajne spremenljivke ne moremo uporabiti analitične inverzne metode. Če računalniški program omogoča numerični izračun inverzne porazdelitvene funkcije normalne porazdelitve, lahko uporabimo inverzno metodo, četudi funkcije F −1(u X i) ne poznamo v eksplicitni obliki. Na primer: v programu EXCEL je funkcija vgrajena kot norminv, v programu MATHEMATICA lahko problem rešimo z ukazom InverseCDF[NormalDistribution[mx,sx]]. Če programska oprema ne omogoča uporabe funkcije F −1(u X i), lahko vzorec normalno porazdeljene slu- čajne spremenljivke generiramo z Box-Müllerjevo polarno metodo. S to metodo vzorec standardizirano normalno porazdeljene slučajne spremenljivke izračunamo iz dveh vzorcev (u1i, i = 1, . . . , n) in (u2i, i = 1, . . . , n) slučajnih spremenljivk U1 in U2, ki sta neodvisni in porazdeljeni enakomerno od 0 do 1, po naslednjih enačbah: p z1i = −2 ln u1i sin(2 π u2i), (7.2) p z2i = −2 ln u1i cos(2 π u2i), kjer z1i in z2i predstavljata vzorec neodvisnih standardizirano normalno porazdeljenih slučajnih spremenljivk Z1 in Z2. Vzorca x1i in x2i poljubnih normalno porazdeljenih slučajnih spremenljivk X1 in X2 s pričakovanima vrednostima mX in m ter standardnima deviacijama σ 1 X2 X1 in σX2 iz vzorcev z1i in z2i določimo s transformacijo: x1i = mX1 + z1i σX1, (7.3) x2i = mX2 + z2i σX2. 92 7 Metoda Monte Carlo 7.1.3 Metoda sprejema/zavrnitve Pri inverzni metodi smo uporabili porazdelitveno funkcijo oziroma inverzno porazdelitveno funkcijo, metoda sprejema in zavrnitve pa temelji na gostoti verjetnosti. Vzemimo, da želimo generirati vzorec slučajne spremenljivke X z zalogo vrednosti na končnem intervalu [a, b]. Taka porazdelitev je na primer porazdelitev beta(4,2)(x − 0.5) za x ∈ [0.5, 1.5], ki je prikazana na sliki 7.3. 2.5 2.0 fX (x) 1.5 1.0 0.5 0.00.0 0.5 1.0 1.5 2.0 x SLIKA 7.3: Gostota verjetnosti, porazdelitev beta(4, 2) Ideja metode sprejema in zavrnitve je v generiranju vzorca enakomerno porazdeljenega slučajnega vektorja na območju x ∈ [a, b] ter y ∈ [0, max (fX (x))], kot je prikazano na sliki 7.4. x∈[a,b] Interval [a, b] je običajno določen z zalogo vrednosti slučajne spremenljivke X. Metoda bi sicer de-lovala tudi, če bi ta interval povečali, a bi s tem povečanjem zmanjšali učinkovitost, saj bi generirali množico točk, za katere zanesljivo vemo, da ne morejo ležati pod grafom gostote verjetnosti in jih zato ne moremo sprejeti. Optimalno je torej najmanjše pravokotno območje, ki očrta graf gostote verjetnosti slučajne spremenljivke X. Točke, ki ležijo pod grafom gostote verjetnosti sprejmemo, to pomeni, da vrednost x dodamo v vzorec. Točke, ki ležijo nad grafom gostote verjetnosti zavrnemo. Na sliki 7.4a sprejete točke označimo z zeleno barvo, zavrnjene pa z rdečo. Na sliki 7.4b vse sprejete točke označimo z zeleno, na tej sliki se vidi, da je rezultat vzorčenja seznam vrednosti na območju od 0.5 do 1.5, pri čemer so točke bolj gosto razporejene tam, kjer je gostota verjetnosti večja. 7.1 Generiranje vzorca slučajne spremenljivke 93 a) b) fX (x) fX (x) x x SLIKA 7.4: Metoda sprejema/zavrnitve Animacija generiranja vzorca slučajne spremenljivke po metodi sprejema in zavrnitve je prikazana na sliki 7.5. 2.5 a) 2.0 fX (x) 1.5 1.0 0.5 0.00.0 0.5 1.0 1.5 2.0 x SLIKA 7.5: Animacija generiranja po metodi sprejema in zavrnitve 7.1.4 Primeri Primer 7.1 Generiranje po inverzni metodi. Z uporabo inverzne metode generirajmo vzorec slučajne spremenljivke T , porazdeljene po eksponentni porazdelitvi, ter slučajne spremenljivke Y , porazdeljene po Gumbelovi porazdelitvi maksimuma. Parameter za eksponentno porazdelitev je: λ = 0.3, za Gumbelovo porazdelitev pa sta: u = 10, α = 0.4. Najprej moramo generirati vzorec enakomerno porazdeljene slučajne spremenljivke U , nato pa z uporabo izrazov v preglednici 7.1 določimo vzorec slučajne spremenljivke T oziroma Y . V preglednici 7.2 prikazujemo prvih deset slučajno generiranih vrednosti ui in iz njih izračunanih vrednosti ti. Na primer: iz prve generirane vrednosti u1 = 0.9384 določimo prvo vrednost v vzorcu eksponentno porazdeljene 94 7 Metoda Monte Carlo slučajne spremenljivke T : 1 1 t1 = − ln(1 − u1) = − ln(1 − 0.9384) = 9.290. λ 0.3 PREGLEDNICA 7.2: Razultat inverze metode – eksponentna porazdelitev Generirane vrednosti enakomerno in eksponentno porazdeljene slučajne spremenljivke ui 0.9384 0.3725 0.5234 0.7788 0.0980 0.8043 0.9363 0.4146 0.8976 0.1660 ti 9.2900 1.5536 2.4700 5.0282 0.3438 5.4369 9.1780 1.7851 7.5975 0.6051 Na sliki 7.6 zgoraj grafično prikazujemo vzorca velikosti n = 50 enakomerno porazdeljene slučajne spremenljivke U in eksponentno porazdeljene slučajne spremenljivke T . Že pri tako majhnem vzorcu lahko prepoznamo lastnosti eksponentne porazdelitve, pri kateri je gostota verjetnosti večja pri nižjih vrednostih spremenljivke in se postopoma znižuje. enakomerna porazdelitev od 0 do 1 eksponentna porazdelitev λ = 0.3 0 2 4 6 8 10 Enakomerna porazdelitev Eksponentna porazdelitev 0.30 1.0 fT(t) fU(u) 0.20 0.6 0.10 0.2 0.0 0.00 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 u t SLIKA 7.6: Inverzna metoda generiranja vzorca – eksponentna porazdelitev Na sliki 7.6 spodaj sta tudi stolpčna diagrama (histograma), ki prikazujeta razporeditev vrednosti večjega vzorca (n = 5000) za obe porazdelitvi. S primerjavo med stolpčnim diagramom in gostoto verjetnosti slučajnih spremenljivk U in T vidimo, da je razporeditev vrednosti ustrezna. V preglednici 7.3 prikazujemo prvih deset slučajno generiranih vrednosti ui in iz njih izračunanih vrednosti yi. Na primer: iz prve generirane vrednosti u1 = 0.4009 določimo prvo vrednost v vzorcu slučajne spremenljivke Y , porazdeljene po Gumbelovi porazdelitvi: 1 1 y1 = u − ln(− ln(u1)) = 10 − ln(− ln(0.4009)) = 10.225. α 0.4 7.1 Generiranje vzorca slučajne spremenljivke 95 PREGLEDNICA 7.3: Razultat inverze metode – Gumbelova porazdelitev Generirane vrednosti enakomerno in Gumbelovo porazdeljene slučajne spremenljivke ui 0.4009 0.7661 0.5612 0.1467 0.9921 0.5106 0.8341 0.7233 0.3591 0.5018 yi 10.225 13.307 11.372 8.3700 22.107 10.993 14.268 12.818 9.9403 10.930 Na sliki 7.7 zgoraj grafično s točkami prikazujemo vzorca velikosti n = 50 enakomerno porazdeljene slučajne spremenljivke U in slučajne spremenljivke Y , porazdeljene po Gumbelovi porazdelitvi. Že pri tako majhnem vzorcu lahko prepoznamo lastnosti Gumbelove porazdelitve, pri kateri je gostota verjetnosti največja v okolici vrednosti u = 10, za manjše in večje vrednosti pa je gostota manjša. Opazimo tudi nesimetričnost porazdelitve. Na isti sliki spodaj sta prikazana tudi stolpčna diagrama za večji vzorec (n = 5000), s katerima ponazorimo razpored vrednosti vzorcev. Le-ta ustrezata izbrani porazdelitvi. enakomerna porazdelitev od 0 do 1 Gumbelova porazdelitev u = 10, α = 0.4 0 5 10 15 20 Enakomerna porazdelitev Gumbelova porazdelitev 0.20 1.0 fY(y) 0.15 fU(u) 0.6 0.10 0.05 0.2 0.0 0.00 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 u y SLIKA 7.7: Inverzna metoda generiranja vzorca – Gumbelova porazdelitev Primer 7.2 Generiranje po metodi Box-Müller. Z uporabo Box-Müllerjeve metode generirajmo vzorec normalno porazdeljene slučajne spremenljivke s srednjo vrednostjo mX = 5 in σX = 2. V preglednici 7.4 prikazujemo prvih deset vrednosti vzorca enakomerno porazdeljenih slučajnih spremenljivk U1 in U2 ter ustreznih vrednosti normalno porazdeljenih slučajnih spremenljivk X1 in X2. Na primer: iz prvih generiranih vrednosti u11 in u21 določimo prvi vrednosti v vzorcu slučajnih spremenljivk X1 in X2: p p x11 = −2 ln(u11) sin(2π u21) = −2 ln(0.3956) sin(2π 0.7620) = 2.2840, (7.4) p p x21 = −2 ln(u11) cos(2π u21) = −2 ln(0.3956) cos(2π 0.7620) = 5.2057. 96 7 Metoda Monte Carlo PREGLEDNICA 7.4: Generiranje po Box-Mülerjevi metodi Generirane vrednosti enakomerno in normalno porazdeljenih slučajnih spremenljivk u1i 0.3956 0.1897 0.7273 0.4939 0.1811 0.5518 0.1928 0.7393 0.4273 0.6994 u2i 0.7620 0.0782 0.1713 0.8727 0.5581 0.4212 0.5647 0.2521 0.000031 0.3338 x1i 2.2840 6.7210 6.4049 3.2959 3.6793 6.0363 3.5654 6.5545 5.0005 6.4624 x2i 5.2047 8.2150 5.7570 6.6549 1.5468 3.0810 1.6667 4.9798 7.6082 4.1503 Na sliki 7.8 s točkami prikazujemo vse štiri vzorce velikosti n = 50 slučajnih spremenljivk U1, U2, X1 in X2. Vidimo, da iz dveh vzorcev slučajnih spremenljivk, porazdaljenih enakomerno od 0 do 1, določimo vzorec dveh neodvisnih spremenljivk, porazdeljenih po normalni porazdelitvi s parametri mX = 5 in σX = 2. S stolpčnima diagramoma prikazujemo še razpored vrednosti v večjih vzorcih n = 5000 in ju primerjamo s teoretičnima porazdelitvama. enakomerna porazdelitev od 0 do 1 normalna porazdelitev m = 5, σ = 2 0 2 4 6 8 10 Enakomerna porazdelitev Normalna porazdelitev 0.25 1.0 0.20 fU1(u) fX1(x) 0.6 0.10 0.2 0.0 0.00 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 u1 x1 SLIKA 7.8: Box-Müllerjeva metoda – normalna porazdelitev Primer 7.3 Uporaba metode sprejema in zavrnitve. Generirajmo vzorec slučajne spremenljivke X, ki je porazdeljena po porazdelitvi, definirani z naslednjim izrazom: 0.275 (x − 1) za 1 ≤ x ≤ 3, fX (x) = 0.1 (6 − x) za 3 < x ≤ 6. Na sliki 7.9 je prikazana gostota verjetnosti slučajne spremenljivke X. 7.1 Generiranje vzorca slučajne spremenljivke 97 0.8 fX (x) 0.6 0.4 0.2 0.00 1 2 3 4 5 6 7 x SLIKA 7.9: Sestavljena trikotna porazdelitev Najprej generiramo dva neodvisna vzorca enakomerno porazdeljenih slučajnih spremenljivk, U1 je enakomerno porazdeljena od 1 do 6, U2 pa od 0 do 0.55. Vzorec 2000 elementov slučajnih spremenljivk U1 in U2 prikazujemo na sliki 7.10. Z zeleno so označene točke, ki jih sprejmemo, s sivo pa tiste, ki jih zavrnemo. 0.8 fX (x) 0.6 0.4 0.2 0.00 1 2 3 4 5 6 7 x SLIKA 7.10: Sestavljena trikotna porazdelitev V preglednici 7.5 je zapisanih prvih desetih generiranih vrednosti slučajnih spremenljivk U1 in U2, izračunane vrednosti gostote verjetnosti ter generirana vrednost vzorca slučajne spremenljivke X, v primeru sprejema in oznaka, če je prišlo do zavrnitve. Na primer: ker je vrednost gostote verjetnosti pri u11 = 2.803 večja od vrednosti u21 = 0.056, vrednost sprejmemo: u21 = 0.056 < fX (u11) = fX (2.803) = 0.496 → x1 = 2.803. Podobno naredimo tudi za drugo in tretjo vrednost. Drugače je pri četrtem elementu vzorca: ker je gostota verjetnosti pri u14 = 3.188 manjša od vrednosti u24 = 0.529, ta element zavrnemo: u24 = 0.529 > fX (u14) = fX (3.453) = 0.255 → ta element zavrnemo. 98 7 Metoda Monte Carlo PREGLEDNICA 7.5: Generiranje po metodi sprejema in zavrnitve Generiranje po metodi sprejema in zavrnitve u1i 2.803 1.850 3.188 3.453 5.587 1.461 3.898 2.983 4.253 1.774 u2i 0.056 0.174 0.040 0.529 0.277 0.160 0.060 0.405 0.449 0.037 fX (u1i) 0.496 0.234 0.281 0.255 0.041 0.127 0.210 0.545 0.175 0.213 xi 2.803 1.850 3.188 - - - 3.898 2.983 - 1.774 Na sliki 7.11 s stolpčnim diagramom prikazujemo razporeditev elementov vzorca in za primerjavo tudi gostoto verjetnosti. Sestavljena trikotna porazdelitev 0.8 fX (x) 0.6 0.4 0.2 0.0 0 1 2 3 4 5 6 7 x SLIKA 7.11: Generirane vrednosti slučajne spremenljivke X 7.2 Generiranje vzorca slučajnega vektorja Pri generiranju vzorca slučajnega vektorja bomo obravnavali tri različne primere: generiranje vzorca poljubno porazdeljenih neodvisnih slučajnih spremenljivk, generiranje vzorca poljubno porazdeljenih odvisnih slučajnih spremenljivk in generiranje vzorca normalno porazdeljenih koreliranih slučajnih spremenljivk. 7.2.1 Neodvisne slučajne spremenljivke – inverzna metoda Generiranje vzorca slučajnega vektorja neodvisnih slučajnih spremenljivk je enako generiranju vzorca posameznih slučajnih spremenljivk. Zato lahko uporabimo poljuben način generiranja, kot sta na primer inverzna metoda ali metoda sprejema in zavrnitve. V primeru uporabe inverzne metode lahko generiranje vzorca povzamemo z naslednjimi izrazi: 7.2 Generiranje vzorca slučajnega vektorja 99 FX (x (u 1 1i) = u1i → x1i = F −1 X 1i), 1 FX (x (u 2 2i) = u2i → x2i = F −1 X 2i), 2 FX (x (u (7.5) 3 3i) = u3i → x3i = F −1 X 3i), 3 · · · FX (x (u k ki) = uki → xki = F −1 X ki), k kjer so x1i, x2i, . . . , xki generirane vrednosti neodvisnih slučajnih spremenljivk X1, X2, . . . , Xk in u1i, u2i, . . . , uki generirane vrednosti neodvisnih slučajnih spremenljivk U1, U2, . . . , Uk, porazdeljene po enakomerni porazdelitvi. Poudariti velja, da moramo pri generiranju vzorca enakomerno porazdeljenih slučajnih spremenljivk Uj le-te generirati kot neodvisne spremenljivke. Za lažje razumevanje generiranja odvisnih slučajnih spremenljivk naredimo še naslednji razmislek, ki bo pripeljal do izrazov (7.5): za neodvisne slučajne spremenljivke velja, da je pogojna porazdelitvena funkcija enaka brezpogojni: P [X1 ≤ x1 ∩ X2 ≤ x2] FX (x1, x2) F 1,X2 X (x = = F (x 2|X1≤x1 2) = P [X2 ≤ x2|X1 ≤ x1] = X 2), P [X 2 1 ≤ x1] FX (x 1 1) (7.6) od koder sledi pogosto uporabljena enačba za določitev porazdelitve slučajnega vektorja neodvisnih spremenljivk, ki pravi, da je porazdelitvena funkcija slučajnega vektorja neodvisnih spremenljivk kar produkt porazdelitvenih funkcij slučajnih spremenljivk: FX (x (x (x 1,X2 1, x2) = FX1 1) FX2 2). (7.7) Ob uporabi inverzne metode generiranja slučajnih spremenljivk porazdelitveno funkcijo slučajnega vektorja {X1, X2} izenačimo s porazdelitveno funkcijo slučajnega vektorja {U1, U2}, ki ga sestavljata neodvisni slučajni spremenljivki, porazdeljeni enakomerno od 0 do 1: FX (x (x (x (u (u (u 1,X2 1, x2) = FX1 1) FX2 2) = FU1,U2 1, u2) = FU1 1) FU2 2) = u1 u2. (7.8) Začnemo z generiranjem vzorca prve slučajne spremenljivke X1: FX (x (u 1 1i) = u1i → x1i = F −1 X 1i). (7.9) 1 Ob upoštevanju (7.9) v (7.8) ostane le izraz, na osnovi katerega lahko generiramo vzorec druge slučajne spremenljivke X2: FX (x (u 2 2i) = u2i → x2i = F −1 X 2i). (7.10) 2 Posplošitev za poljubno število neodvisnih slučajnih spremenljivk je preprosta. 100 7 Metoda Monte Carlo 7.2.2 Odvisne slučajne spremenljivke – inverzna metoda Za dve odvisni slučajni spremenljivki X1 in X2 velja, da porazdelitvene funkcije slučajnega vektorja {X1, X2} ne moremo zapisati kot produkt porazdelitvenih funkcij FX (x (x 1 1) in FX2 2), kot pri neodvis- nih slučajnih spremenljivkah. Velja namreč: P [X1 ≤ x1 ∩ X2 ≤ x2] FX (x1, x2) F 1,X2 X (x = , (7.11) 2|X1≤x1 2) = P [X2 ≤ x2|X1 ≤ x1] = P [X1 ≤ x1] FX (x 1 1) od koder sledi zveza: FX (x (x (x (x (x 1,X2 1, x2) = FX1 1) FX2|X1≤x1 2) ̸= FX1 1) FX2 2). (7.12) Zato moramo vzorec slučajnih spremenljivk generirati postopoma. Najprej generiramo vzorec prve slučajne spremenljivke: x1i = F −1(u X 1i). (7.13) 1 Nato poiščemo pogojno porazdelitveno funkcijo FX (x 2|X1=x1 2), na osnovi katere bomo lahko generi- rali vzorec naslednje slučajne spremenljivke. To porazdelitveno funkcijo najlažje določimo iz pogojne gostote verjetnosti, ki jo določimo po naslednji enačbi: fX (x1, x2) f 1,X2 X (x . (7.14) 2|X1=x1 2) = fX (x 1 1) Porazdelitveno funkcijo določimo z integracijo gostote verjetnosti: Z x2 FX (x f (˜ x) d˜ x. (7.15) 2|X1=x1 2) = X2|X1=x1 −∞ Vzorec slučajne spremenljivke x2 generiramo z uporabo inverzne porazdelitvene funkcije (7.15): FX (x (u 2|X1=x1i 2i) = u2i → x2i = F −1 X 2i). (7.16) 2|X1=x1i Podobno naredimo tudi za vse nadaljnje spremenljivke: fX (x1, x2, x3) f 1,X2,X3 X (x , (7.17) 3|X1=x1∩X2=x2 3) = fX (x 1,X2 1, x2) opravimo integracijo za določitev pogojne porazdelitvene funkcije: Z x3 FX (x f (˜ x) d˜ x (7.18) 3|X1=x1∩X2=x2 3) = X3|X1=x1∩X2=x2 −∞ in vzorec slučajne spremenljivke x3 generiramo z uporabo inverzne porazdelitvene funkcije (7.18): FX (x (u 3|X1=x1i∩X2=x2i 3i) = u3i → x3i = F −1 X 3i). (7.19) 3|X1=x1i∩X2=x2i 7.2 Generiranje vzorca slučajnega vektorja 101 Na osnovi enačb (7.13), (7.16) in (7.19) lahko zapišemo izraze za generiranje vzorca poljubno porazdeljenih odvisnih slučajnih spremenljivk: x1i = F −1(u X 1i), 1 x2i = F −1 (u X 2i), 2|X1=x1i x3i = F −1 (u (7.20) X 3i), 3|X1=x1i∩X2=x2i · · · xki = F −1 (u X ki), k |X1=x1i∩X2=x2i∩...∩Xk−1=xk−1,i kjer so x1i, x2i, . . . , xki generirane vrednosti neodvisnih slučajnih spremenljivk X1, X2, . . . , Xk, in u1i, u2i, . . . , uki generirane vrednosti neodvisnih slučajnih spremenljivk U1, U2, . . . , Uk, porazdeljene po enakomerni porazdelitvi. Če primerjamo enačbe (7.5) in (7.20), opazimo, da za generiranje vzorca odvisne slučajne spremenljivke X2 potrebujemo vzorčno vrednost slučajne spremenljivke x1i, za generiranje vzorca X3 potrebujemo vzorčne vrednosti x1i in x2i in tako dalje, medtem ko pri generiranju vzorca neodvisnih slučajnih spremenljivk teh vrednosti ne potrebujemo. 7.2.3 Korelirane normalno porazdeljene spremenljivke Pri generiranju slučajnega vektorja X koreliranih normalno porazdeljenih slučajnih spremenljivk bomo izhajali iz vektorja Y standardizirano normalnih neodvisnih slučajnih spremenljivk. Vektor X je določen z vektorjem srednjih vrednosti mX in kovariančno matriko CX, medtem ko so srednje vrednosti slučajnega vektorja Y enake nič in je kovariančna matrika enotska matrika CY = I:  m    X σ2 σ · · · σ 1 X X X 1 1X2 1Xn  mX σ2 σ2 · · · σ 2   X X X2Xn  m 2X1 2     X = . , CX = . . . , (7.21)  . . .   . . . . .     . . .  mX σ σ · · · σ2 n XnX1 XnX2 Xn  0   1 0 · · · 0   0   0 1 · · · 0  m     Y = . , CY = I = . . . . . (7.22)  ..   .. .. . . ..      0 0 0 · · · 1 Med njima velja linearna zveza, ki jo v matrični obliki zapišemo z naslednjim izrazom: X = mX + HY → X − mX = HY, (7.23) kjer je H trikotna transformacijska matrika, ki jo moramo še določiti. Za določitev matrike H izhajamo iz enačbe za kovariančno matriko, v kateri upoštevamo, da je varianca oziroma kovarianca pričakovana vrednost funkcije slučajne spremenljivke: σ2X = E[(X1 − m )2] oziroma σ = E[(X1 − m )(X2 − m )]. (7.24) 1 X1 X1X2 X1 X2 102 7 Metoda Monte Carlo V matrični obliki kovariančno matriko zapišemo takole: h CX = E (X − mX)(X − mX)T i . (7.25) Če v enačbo (7.25) vstavimo (7.23), dobimo: h h i CX = E (HY)(HY)T i = E H(YYT )HT = ... upoštevamo, da matrika H ni slučajna, zato jo lahko izpostavimo ... = H E YYT HT = H C (7.26) Y HT = ... upoštevamo, da je kovariančna matrika enotska: CY = I ... = H I HT = H HT . V primeru, da simetrično pozitivno definitno matriko zapišemo kot produkt dveh matrik H in HT , lahko matriko H določimo z razcepom Choleskega1. 7.2.4 Primeri Primer 7.4 Generiranje po inverzni metodi. Z uporabo inverzne metode generirajmo vzorec slučajnega vektorja {X, Y }, porazdeljenega po porazdelitvi, definirani z naslednjim izrazom: fXY (x, y) = 3 x, za 0 ≤ x ≤ 1, x − 1 ≤ y ≤ 1 − x. (7.27) Grafa gostote verjetnosti in porazdelitvene funkcije sta prikazana na sliki 7.12. Gostota verjetnosti Porazdelitvena funkcija 1 0 -1 1 SLIKA 7.12: Gostota verjetnosti in porazdelitvena funkcija slučajnega vektorja {X, Y } Določimo najprej robni gostoti verjetnosti fX (x) in fY (y): Z ∞ Z 1−x 1−x f X (x) = fXY (x, y) dy = 3 x dy = 3 x y = 6x(1 − x), (7.28) −∞ x−1 x−1 1 André-Louis Cholesky, francoski oficir in matematik (1875–1918). 7.2 Generiranje vzorca slučajnega vektorja 103  Z 1+y 3x2 1+y 3(1 + y)2  3 x dx = = za y ≤ 0, Z ∞   2 0 2 f 0 Y (y) = fXY (x, y) dx = (7.29) Z 1−y −∞ 3x2 1−y 3(1 − y)2    3 x dx = = za y > 0. 0 2 0 2 Grafa obeh robnih gostot verjetnosti prikazujemo na sliki 7.13. Za generiranje po inverzni metodi potre-1.5 1.5 f fY (y) X (x) 1.0 1.0 0.5 0.5 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 - 1.0 - 0.5 0.0 0.5 1.0 x y SLIKA 7.13: Robni gostoti verjetnosti fX (x) in fY (y) bujemo inverzne porazdelitvene funkcije, kot je zapisano z enačbami (7.20). Najprej določimo inverzno porazdelitveno funkcijo slučajne spremenljivke X: Z x Z x FX (x) = fX (˜ x) d˜ x = 6˜ x(1 − ˜ x) d˜ x = 3x2 − 2x3, −∞ 0 ... pri določitvi inverzne funkcije kubičnega polinoma smo upoštevali le eno rešitev, tisto, ki ustreza vrednostim x ∈ [0, 1] v območju zaloge vrednosti X ... (7.30)  √ 2/3 √  1 1 + i 3 2u + 2p(u − 1)u − 1 + 1 − i 3 F −1(u) =  + 2 . X 4  1/3  2u + 2p(u − 1)u − 1 Na prvi pogled zadnji izraz predstavlja kompleksno število, a se izkaže, da je za vrednosti u ∈ [0, 1] rezultat realno število med 0 in 1. V primeru, da analitične rešitve inverzne funkcije ne moremo najti, lahko uporabimo numerične metode. Sedaj moramo določiti še pogojno gostoto verjetnosti fY |X=x(y) po naslednji enačbi; f 3x XY (x, y) = 1 za x − 1 ≤ y ≤ 1 − x, f 6x(1−x) 2(1−x) Y |X=x(y) = = (7.31) fX (x) 0 drugje. Ugotovimo lahko, da je pogojna slučajna spremenljivka Y pri določeni vrednosti x enakomerno porazdeljena od x − 1 do 1 − x. Določitev te pogojne porazdelitve grafično prikazujemo na sliki 7.14. 104 7 Metoda Monte Carlo z 2.0 x = 0.7 1 x = 0.3 1.5 0 x = 0.7 fY X=x( y) 1.0 -1 x = 0.3 0.5 0.0 - 1.0 - 0.5 0.0 0.5 1.0 y SLIKA 7.14: Določitev pogojne gostote verjetnosti fY |X=x(y) za dve vrednosti x (0.3 in 0.7) Če želimo uporabiti inverzno metodo generiranja, moramo določiti še pogojno porazdelitveno funkcijo in njeno inverzno funkcijo:  0 za y < x − 1, Z y  F y−x+1 Y |X=x(y) = fY |X=x(˜ y) d˜ y = za x − 1 ≤ y ≤ 1 − x, 2(1−x) −∞  1 za y > 1 − x, (7.32) ... poiščemo inverzno porazdelitveno funkcijo: F −1 (u) = 2(1 − x)u − (1 − x). Y |X=x Pripravimo enačbi inverznih porazdelitvenih funkcij tako, da bomo lahko uporabili izraze (7.20) za generiranje odvisnih slučajnih spremenljivk:  √ 2/3 √  1 1 + i 3 2u1i + 2p(u1i − 1)u1i − 1 + 1 − i 3 x   i = + 2 , 4  1/3  (7.33) 2u1i + 2p(u1i − 1)u1i − 1 yi = 2(1 − xi)u2i − (1 − xi). V preglednici 7.6 prikazujemo vrednosti prvih desetih elementov vzorca. Na primer, iz prve vrednosti u11 = 0.076 po prvi od enačb (7.33) izračunamo prvo vrednost x1 = 0.169. Nato uporabimo drugo izmed enačbo (7.33) in iz u21 = 0.561 in x1 = 0.169 izračunamo še y1 = 0.101. PREGLEDNICA 7.6: Generiranje slučajnega vektorja po inverzni metodi Generirane vrednosti u1i 0.076 0.806 0.101 0.405 0.202 0.595 0.709 0.289 0.185 0.871 u2i 0.561 0.201 0.732 0.221 0.240 0.894 0.880 0.658 0.914 0.062 xi 0.169 0.718 0.196 0.436 0.289 0.564 0.644 0.355 0.275 0.775 yi 0.101 −0.169 0.373 −0.315 −0.370 0.344 0.271 0.203 0.601 −0.197 7.2 Generiranje vzorca slučajnega vektorja 105 Na sliki 7.15 prikazujemo vzorca 10000 elementov slučajnih vektorjev {U1, U2} in {X, Y }. Vidimo razliko med enakomerno porazdeljenim vektorjem {U1, U2} in neenakomerno porazdelitvijo vektorja {X, Y }, ki se odraža v povečani gostoti točk pri višjih vrednostih X. Enakomerna porazdelitev 1.2 Slučajni vektor {X,Y } 1.0 1.0 y u2 0.5 0.8 0.6 0.0 0.4 - 0.5 0.2 - 1.0 0.00.0 0.2 0.4 0.6 0.8 1.0 1.2 u 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1 x SLIKA 7.15: Vzorec slučajnih vektorjev {U1, U2} in {X, Y } Na naslednji sliki 7.16 prikazujemo stolpčne diagrame vzorcev slučajnih spremenljivk X in Y ter ju primerjamo s teoretičnima gostotama porazdelitve. Očitno je, da smo obe slučajni spremenljivki generirali skladno z zahtevano porazdelitvijo. Slučajna spremenljivka Slučajna spremenljivka X Y 1.5 1.5 f f Y ( y) X (x) 1.0 1.0 0.5 0.5 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 - 1.0 - 0.5 0.0 0.5 1.0 x y SLIKA 7.16: Vzorec slučajnih spremenljivk in primerjava z njunima gostotama verjetnosti Na zadnji sliki 7.17 prikazujemo še stolpčni diagram za vzorec slučajnega vektorja {X, Y } in ga primerjamo s teoretično gostoto verjetnosti. Tudi tukaj lahko ugotovimo, da se generirani vzorec sklada s predpostavljeno porazdelitvijo. 106 7 Metoda Monte Carlo 1 1 -1 -1 0 0 1 1 SLIKA 7.17: Vzorec slučajnega vektorja {X, Y } in primerjava z njegovo gostoto verjetnosti Primer 7.5 Generiranje po metodi sprejema in zavrnitve. Z uporabo metode sprejema in zavrnitve generirajmo vzorec slučajnega vektorja {X, Y } iz primera 7.4. Pri uporabi te metode moramo generirati tri neodvisne enakomerno porazdeljene slučajne spremenljivke. Glede na zalogo vrednosti slučajnega vektorja {X, Y } je prva slučajna spremenljivka U1 enakomerno porazdeljena od 0 do 1, druga U2 je enakomerno porazdeljena od −1 do 1, glede na največjo vrednost gostote verjetnosti pa je tretja slučajna spremenljivka U3 porazdeljena enakomerno od 0 do 3. Na sliki 7.18 prikazujemo množico enakomerno porazdeljenih točk v kvadru. Zeleno so obarvane tiste točke, ki jih sprejmemo. Rezultat generiranja po metodi sprejema in zavrnitve je prikazan na desni strani iste slike. V preglednici 7.7 prikazujemo določitev vzorca slučajnega vektorja po metodi sprejema in zavrnitve. Prikazujemo prvih deset slučajno izbranih točk. Za prvo točko {0.367, −0.161, 0.736} lahko ugotovimo, da je gostota verjetnosti: u31 = 0.736 < fXY (u11, u21) = fXY (0.367, −0.161) = 1.101, zato točko sprejmemo. V vzorec dodamo par {0.367, −0.161}. Naslednjo točko {0.891, 0.384, 2.640} moramo žal zavrniti, saj je za koordinati {0.891, 0.384} gostota verjetnosti enaka nič (glej enačbo (7.27)): u32 = 2.640 > fXY (u12, u22) = fXY (0.891, 0.384) = 0. Vidimo, da lahko od prvih desetih generiranih točk sprejmemo le dve (prvo in šesto), druge pa moramo zavrniti. Delež točk, ki jih bomo sprejeli, običajno lahko ocenimo s primerjavo prostornin: prostornina kvadra, v katerem generiramo točke, je 6, prostornina piramide, ki ustreza območju sprejema, pa je enaka 1. Zato lahko pričakujemo, da bo delež sprejetih točk okoli 1/6. Za vzorec 10000 elementov smo v naši simulaciji generirali 60066 točk, kar potrjuje našo oceno. 7.2 Generiranje vzorca slučajnega vektorja 107 PREGLEDNICA 7.7: Generiranje slučajnega vektorja po metodi sprejema in zavrnitve Generirane vrednosti u1i 0.367 0.891 0.527 0.678 0.397 0.542 0.603 0.894 0.353 0.370 u2i −0.161 0.384 −0.800 −0.575 0.763 −0.029 0.562 0.729 −0.051 −0.910 u3i 0.736 2.640 0.563 1.830 2.204 1.457 2.158 0.784 2.670 1.783 fXY 1.101 0 0 0 0 1.627 0 0 1.058 0 xi 0.367 - - - - 0.542 - - - - yi −0.161 - - - - −0.029 - - - - Slučajni vektor {X,Y } 1.0 y 0.5 0.0 - 0.5 - 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 x SLIKA 7.18: Generiranje množice enakomerno razporejenih točk po kvadru Primer 7.6 Generiranje koreliranih normalno porazdeljenih slučajnih vektorjev. Z uporabo opisane metode in razcepa Choleskega generirajmo slučajni vektor X z naslednjimi lastnostmi: m 10 σ2 σ 12 6 m X1 X X 1 1X2 X = = , C = . (7.34) m X = X 15 σ σ2 6 6 2 X1X2 X2 Gostoto verjetnosti slučajnega vektorja X in podobnega vektorja z neodvisnimi slučajnimi spremenljivkami (σX = 0) prikazujemo na sliki 7.19. 1X2 108 7 Metoda Monte Carlo Korelirani vektor x Nekorelirani vektor 2 x 2 fX X f 1 2 X X 1 2 x 1 x 1 SLIKA 7.19: Gostota verjetnosti vektorja X koreliranih in nekoreliranih slučajnih spremenljivk Pred uporabo enačbe (7.23) moramo določiti matriko H, ki jo dobimo z razcepom Choleskega matrike CX. Matrika H je: √ 2 3 0 3.464 0 H = √ √ = . (7.35) 3 3 1.732 1.732 Nato generiramo vzorec neodvisnih standardizirano normalno porazdeljenih slučajnih spremenljivk Y1 in Y2. Prvih deset slučajno generiranih vrednosti slučajnih vektorjev Y in X prikazujemo v preglednici 7.8. Na primer, za prvi element vzorca izračunamo: 10 3.464 0 0.548 11.90 X1 = + = . (7.36) 15 1.732 1.732 1.287 18.19 Na podoben način izračunamo tudi vse druge elemente vzorca slučajnega vektorja X. PREGLEDNICA 7.8: Generiranje odvisnih normalno porazdeljenih slučajnih spremenljivk Generirane vrednosti y1i 0.548 −0.416 −1.664 0.392 −1.857 1.984 0.912 1.278 −1.019 −0.738 y2i 1.287 0.653 −1.178 1.575 0.763 −1.082 0.751 0.738 1.813 −0.245 x1i 11.90 8.56 4.23 11.36 3.57 16.87 13.16 14.43 6.47 7.44 x2i 18.19 15.41 10.08 18.41 13.11 16.56 17.88 18.49 16.38 13.30 Na sliki 7.20 prikazujemo vzorca slučajnega vektorja neodvisnih oziroma nekoreliranih standardizirano normalno porazdeljenih spremenljivk Yi ter odvisnih oziroma koreliranih normalno porazdeljenih slu- čajnih spremenljivk Xi. 7.2 Generiranje vzorca slučajnega vektorja 109 Neodvisni standardizirano normalni Odvisni normalni slučajni spremenljivki 3 25 y2 x2 2 20 1 15 0 10 - 1 5 - 2 - 3 0 - 3 - 2 - 1 0 1 2 3 y 0 5 10 15 20 25 1 x1 SLIKA 7.20: Vzorca slučajnih vektorjev Y in X Če bi želeli generirati neodvisne normalno porazdeljene slučajne spremenljivke, bi uporabili enačbo: xki = mX + σ y k Xk ki. (7.37) Prvih deset naključno izbranih elementov vzorca prikazujemo v preglednici 7.9. S primerjavo med pre-glednicama 7.8 in 7.9 vidimo, da se vrednosti vzorca slučajne spremenljivke X1 ne razlikujejo, do razlik pride pri generiranju druge slučajne spremenljivke X2, ki je v zadnjem primeru neodvisna od slučajne spremenljivke X1. PREGLEDNICA 7.9: Generiranje neodvisnih normalno porazdeljenih slučajnih spremenljivk Generirane vrednosti y1i 0.548 −0.416 −1.664 0.392 −1.857 1.984 0.912 1.278 −1.019 −0.738 y2i 1.287 0.653 −1.178 1.575 0.763 −1.082 0.751 0.738 1.813 −0.245 x1i 11.90 8.56 4.23 11.36 3.57 16.87 13.16 14.43 6.47 7.44 x2i 18.15 16.60 12.11 18.86 16.87 12.35 16.84 16.81 19.44 14.40 Na sliki 7.21 prikazujemo vzorca slučajnega vektorja neodvisnih oziroma nekoreliranih standardizirano normalno porazdeljenih spremenljivk Yi ter neodvisnih poljubno normalno porazdeljenih slučajnih spremenljivk Xi. V primerjavi med slikama 7.20 in 7.21 vidimo, da sta pri koreliranih slučajnih spremenljivkah osi oblaka točk poševni, medtem ko sta pri nekoreliranih osi oblaka točk vzporedni s koordinatnima osema. 110 7 Metoda Monte Carlo Neodvisni standardizirano normalni Neodvisni normalni slučajni spremenljivki 3 25 y2 x2 2 20 1 15 0 10 - 1 5 - 2 - 3 0 - 3 - 2 - 1 0 1 2 3 y 0 5 10 15 20 25 1 x1 SLIKA 7.21: Vzorca slučajnih vektorjev Y in X Primer 7.7 Generiranje koreliranih normalno porazdeljenih slučajnih vektorjev. Z uporabo opisane metode in razcepa Choleskega generirajmo slučajni vektor štirih slučajnih spremenljivk X z naslednjimi lastnostmi:  m        X 10 σ2 σ σ σ 12 6 4 2 1 X X X X 1 1X2 1X3 1X4 m 15 σ σ2 σ σ 6 6 4 4 m  X2     X1X2 X X X    2 2X3 2X4 X =   =   , CX =   =   .  mX 10 σ σ σ2 σ 4 4 10 4 3     X1X3 X2X3 X X    3 3X4 mX 20 σ σ σ σ2 2 4 4 8 4 X1X4 X2X4 X3X4 X4 Najprej moramo določiti matriko H, ki jo dobimo z razcepom Choleskega matrike CX: √  2 3 0 0 0  √ √  3.464 0 0 0   3 3 0 0  1.732 1.732 0 0 H =  q     2 = . √ 2 √ 22 0    (7.38)  1.155 1.155 2.708 0 3 3 3     √ q q  1 √ 3 2 2 146 0.577 1.732 0.492 2.103 3 33 33 V preglednici 7.10 prikazujemo prvih osem naključno izbranih elementov vzorcev slučajnih vektorjev Y in X. Na primer, tretji element vzorca slučajnega vektorja X določimo takole:  10   3.464 0 0 0   0.723   12.50  15 1.732 1.732 0 0 0.200 16.60 X         3 =   +    = . (7.39) −     10   1.155 1.155 2.708 0   1.182   7.87  20 0.577 1.732 0.492 2.103 −1.784 16.43 7.2 Generiranje vzorca slučajnega vektorja 111 PREGLEDNICA 7.10: Generiranje odvisnih normalno porazdeljenih slučajnih spremenljivk Generirane vrednosti y1i −0.644 0.555 0.723 −1.389 0.762 −1.429 1.610 0.627 y2i 0.286 0.765 0.200 −0.009 −0.244 −0.101 −1.416 0.739 y3i 1.582 −0.668 −1.182 1.026 −0.594 −1.274 −1.671 2.605 y4i 0.069 −0.960 −1.784 −1.069 −0.532 −0.938 −0.914 −0.270 x1i 7.77 11.92 12.50 5.19 12.64 5.05 15.58 12.17 x2i 14.38 17.29 16.60 12.58 15.90 12.35 15.34 17.37 x3i 13.87 9.72 7.87 11.16 8.99 4.78 5.70 18.63 x4i 21.05 19.30 16.43 17.44 18.61 16.40 15.73 22.36 Na sliki 7.22 prikazujemo vzorec slučajnega vektorja X. Ker je vektor štirirazsežen, grafično prikazujemo posamezne dele vektorja in sicer po parih slučajnih spremenljivk. Iz slike lahko vidimo, da je korelacija med X1 in X2 močnejša kot korelacija med X1 in X3 ali X1 in X4. To se seveda odraža že iz podatkov, saj je korelacijski koeficient ρX = 0.707 > ρ = 0.365 > ρ = 0.204. 1X2 X1X3 X1X4 30 30 30 x2 25 x3 25 x4 25 20 20 20 15 15 15 10 10 10 5 5 5 0 0 0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 x1 x1 x1 30 30 30 x 3 25 x x 4 25 4 25 20 20 20 15 15 15 10 10 10 5 5 5 0 0 0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 x2 x2 x3 SLIKA 7.22: Vzorec slučajnega vektorja X 112 7 Metoda Monte Carlo 8 Literatura 1. R. Andričević, H. Gotovac, I. Ljubenkov, Geostatistika, umijeće prostorne analize, Gradževinsko-Arhitektosnki Fakultet Sveučilišta u Splitu, Split, 2007. 2. J.-P. Chilès, P. Delfiner, Geostatistics, Modeling Spatial Uncertainty, John Wiley & Sons, New York, 1999. 3. N. Christou, Universal kriging, Class notes, University of California, Los Angeles, 2014. 4. N. E. Cressie, Statistics for Spatial Data, John Wiley & Sons, New York, 1993. 5. T. Hengl, A Practical Guide to Geostatistical Mapping, Office for Official Publications of the European Communities, Luxembourg, 2009. 6. K, Krivoruchko, Spatial Statistical Data Analysis for GIS Users, Esri Press, Redlands, California, 2011. 7. A. Paez, D. M. Scott, Spatial statistics for urban analysis: A review of techniques with examples, GeoJournal, let. 61, str. 53–67, 2004. 8. R. Webster, O. Atteia, J.P. Dubois, Coregionalization of trace metals in the soil in the Swiss Jura, European Journal of Soil Science, let. 45, str. 205–218, 1994. Stvarno kazalo aproksimacija variograma, 31 običajno, 50, 57 preprosto, 46, 56 Box–Müllerjeva metoda, 91 s trendom, 51, 58 krigiranje s trendom, 51 eksponentna porazdelitev, 90 kurtosis, 12 enakomerna porazdelitev, 89, 90 Excel, 91 Lagrangev množitelj, 50, 52 linearna regresija, 61 Fréchetova porazdelitev, 90 Mathematica, 91 Gearyjevo razmerje, 36, 40 metoda kvadratov, 70 generiranje vzorca slučajne spremenljivke, 89 metoda Monte Carlo, 87 Box–Müllerjeva metoda, 91 metoda najbližjega soseda, 70, 77 inverzna metoda, 89, 98, 100 metoda najmanjših kvadratov, 62, 63 metoda sprejema/zavrnitve, 92, 93 momenti slučajnih spremenljivk, 12 normalna porazdelitev, 91, 101 Moranov indeks, 36, 37 generiranje vzorca slučajnega vektorja, 98 geostatistični podatki, 4, 23 neodvisnost, 11 Jura, 5, 24 normalna porazdelitev, 91, 107 gostota verjetnosti, 10 Gumbelova porazdelitev, 90 običajno krigiranje, 50 odstopanje, 61 indikatorji, 33 osamelec, 62 koeficient simetričnosti, 12 podatki v mreži, 6 korelacijska funkcija, 14 Poissonova porazdelitev, 71 korelacijski koeficient, 13 porazdelitev kovariančna funkcija, 14, 15 eksponentna, 90 ocena, 29 enakomerna, 89, 90 kovariančna funkcija z indikatorji, 34 Fréchetova, 90 kovarianca, 13 Gumbelova, 90 kovariogram, 20 normalna, 91, 107 križna kovariančna funkcija, 20 Weibullova, 90 krigiranje, 45 porazdelitvena funkcija, 9 115 116 Stvarno kazalo pogojna, 11 prag, 17 robna, 11 učinek zrna, 17 preizkušanje domnev, 35 variogram z indikatorji, 34 preizkus χ2, 72 vektorsko slučajno polje preprosto krigiranje, 46 kovariogram, 20 pričakovana vrednost, 12 križna kovariančna funkcija, 20 prostorski vzorci, 8, 70 križni variogram, 20 prostorsko utežena regresija, 61 porazdelitvena funkcija, 19 verjetnostna funkcija, 10 razcep Choleskega, 102, 107 vezani ekstrem, 50, 52 razsevni graf, 28 vzorčna kovariančna funkcija, 29 regresija vzorčni variogram, 30 linearna, 61 rezidual, 61 Weibullova porazdelitev, 90 robustna metoda, 62 zaloga vrednosti, 9 semivariogram, 15 slučajna funkcija, 13 slučajna spremenljivka, 9 slučajni proces, 13 slučajni vektor, 10 slučajno polje, 13 gostota verjetnosti, 14, 19 korelacijska funkcija, 14 kovariančna funkcija, 14 porazdelitvena funkcija, 13 srednja vrednost, 14, 19 varianca, 14, 19 variogram, 15 srednja vrednost slučajnega polja, 14, 19 stacionarno polje kovariančna funkcija, 15 variogram, 15 stacionarnost, 15 testiranje hipotez, 35 varianca, 12 varianca slučajnega polja, 14, 19 variogram, 15 doseg, 17 empirična funkcija, 31 ocena, 30 Document Outline Predgovor Vsebina Uvod Zgodovina prostorske statistike Podatki za prostorsko statistiko Vrste obravnavnih slučajnih spremenljivk Geostatistični podatki Podatki v mreži Prostorski točkovni vzorci Slučajne funkcije, polja in procesi Slučajne spremenljivke Slučajni vektorji Momenti slučajnih spremenljivk in vektorjev Slučajne funkcije Vektorske slučajne funkcije Geostatistika Geostatistični podatki Priprava podatkov Razsevni graf Kovariančna funkcija Variogram Določitev empirične funkcije variograma Geostatistična analiza z indikatorji Preizkušanje domnev o prostorski povezanosti Moranov indeks Gearyjevo razmerje Krigiranje Preprosto krigiranje Običajno krigiranje Krigiranje s trendom Primeri krigiranja Preprosto krigiranje Običajno krigiranje Krigiranje s trendom Prostorsko utežena regresija Linearna regresija Utežena linearna regresija Določitev uteži Izbira funkcije za določitev uteži in kritične razdalje Primeri Prostorski vzorci Metoda kvadratov Primeri – metoda kvadratov Pomanjkljivosti metode kvadratov Metoda najbližjega soseda Primeri – metoda najbližjega soseda Metoda Monte Carlo Generiranje vzorca slučajne spremenljivke Inverzna metoda generiranja vzorca Generiranje vzorca normalno porazdeljenih slučajnih spremenljivk Metoda sprejema/zavrnitve Primeri Generiranje vzorca slučajnega vektorja Neodvisne slučajne spremenljivke – inverzna metoda Odvisne slučajne spremenljivke – inverzna metoda Korelirane normalno porazdeljene spremenljivke Primeri Literatura Stvarno kazalo Prazna stran Prazna stran