13 Symulacje

Odtworzenie obliczeń z tego rozdziału wymaga załączenia poniższych pakietów oraz wczytania poniższych danych:

13.1 Symulacje geostatystyczne

Kriging daje optymalne estymacje, czyli wyznacza najbardziej potencjalnie możliwą wartość dla wybranej lokalizacji. Dodatkowo, efektem krigingu jest wygładzony obraz. W konsekwencji wyniki estymacji różnią się od danych pomiarowych. Uzyskiwana jest tylko (aż?) estymacja, a prawdziwa wartość jest niepewna… Korzystając z symulacji geostatystycznych nie tworzymy estymacji, ale generujemy równie prawdopodobne możliwości poprzez symulację z rozkładu prawdopodobieństwa (wykorzystując generator liczb losowych).

13.1.1 Właściwości

Właściwości symulacji geostatystycznych:

  • Efekt symulacji ma bardziej realistyczny przestrzenny wzór (ang. pattern) niż kriging, którego efektem jest wygładzona reprezentacja rzeczywistości
  • Każda z symulowanych map jest równie prawdopodobna
  • Symulacje pozwalają na przedstawianie niepewności interpolacji
  • Jednocześnie - kriging jest znacznie lepszy, gdy naszym celem jest jak najdokładniejsza estymacja

13.1.2 Typy symulacji

Istnieją dwa typy symulacji geostatystycznych:

  • Symulacje bezwarunkowe (ang. Unconditional Simulations) - wykorzystujące semiwariogram, żeby włączyć informację przestrzenną, ale wartości ze zmierzonych punktów nie są w niej wykorzystywane
  • Symulacje warunkowe (ang. Conditional Simulations) - opiera się ona o średnią wartość, strukturę kowariancji oraz obserwowane wartości (w tym sekwencyjna symulacja danych kodowanych (ang. Sequential indicator simulation))

13.2 Symulacje bezwarunkowe

Symulacje bezwarunkowe w pakiecie gstat tworzy się z wykorzystaniem funkcji krige(). Podobnie jak w przypadku estymacji geostatystycznych, należy tutaj podać wzór, model, siatkę, średnią globalną (beta), oraz liczbę sąsiednich punktów używanych do symulacji (w przykładzie poniżej nmax = 30). Należy wprowadzić również informacje, że nie korzystamy z danych punktowych (locations = NULL) oraz że chcemy stworzyć dane sztuczne (dummy = TRUE). Ostatni argument (nsim = 4) informuje o liczbie symulacji do przeprowadzenia.

## [using unconditional Gaussian simulation]

## [using unconditional Gaussian simulation]

13.3 Symulacje warunkowe

Jednym z podstawowych typów symulacji warunkowych jest sekwencyjna symulacja gaussowska (ang. Sequential Gaussian simulation). Polega ona na:

  1. Wybraniu lokalizacji nie posiadającej zmierzonej wartości badanej zmiennej
  2. Krigingu wartości tej lokalizacji korzystając z dostępnych danych, co pozwala na uzyskanie rozkładu prawdopodobieństwa badanej zmiennej
  3. Wylosowaniu wartości z rozkładu prawdopodobieństwa za pomocą generatora liczba losowych i przypisaniu tej wartości do lokalizacji
  4. Dodaniu symulowanej wartości do zbioru danych i przejściu do kolejnej lokalizacji
  5. Powtórzeniu poprzednich kroków, aż do momentu gdy nie pozostanie już żadna nieokreślona lokalizacja

Sekwencyjna symulacja gaussowska wymaga zmiennej posiadającej wartości o rozkładzie zbliżonym do normalnego. Można to sprawdzić poprzez wizualizacje danych (histogram, wykres kwantyl-kwantyl) lub też test statystyczny (test Shapiro-Wilka). Zmienna temp nie ma rozkładu zbliżonego do normalnego.

## 
##  Shapiro-Wilk normality test
## 
## data:  punkty$temp
## W = 0.97683, p-value = 0.0004194

Na potrzeby symulacji zmienna temp została zlogarytmizowna.

## 
##  Shapiro-Wilk normality test
## 
## data:  punkty$temp_log
## W = 0.98907, p-value = 0.05568

Dalsze etapy przypominają przeprowadzenie estymacji statystycznej, jedynym wyjątkiem jest dodanie argumentu mówiącego o liczbie symulacji do przeprowadzenia (nsim w funkcji krige()).

## drawing 4 GLS realisations of beta...
## [using conditional Gaussian simulation]

Wyniki symulacji można przetworzyć do pierwotnej jednostki z użyciem funkcji wykładniczej (exp).

## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min      max
## x 745541.7 756971.7
## y 712616.2 721256.2
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Number of points: 10993
## Grid attributes:
##    cellcentre.offset cellsize cells.dim
## s1          745586.7       90       127
## s2          712661.2       90        96
## Data attributes:
##       sim1            sim2            sim3            sim4      
##  Min.   :2.058   Min.   :2.009   Min.   :2.053   Min.   :1.994  
##  1st Qu.:2.560   1st Qu.:2.559   1st Qu.:2.556   1st Qu.:2.547  
##  Median :2.725   Median :2.717   Median :2.714   Median :2.720  
##  Mean   :2.733   Mean   :2.735   Mean   :2.724   Mean   :2.726  
##  3rd Qu.:2.902   3rd Qu.:2.919   3rd Qu.:2.889   3rd Qu.:2.903  
##  Max.   :3.473   Max.   :3.434   Max.   :3.376   Max.   :3.465
## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min      max
## x 745541.7 756971.7
## y 712616.2 721256.2
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Number of points: 10993
## Grid attributes:
##    cellcentre.offset cellsize cells.dim
## s1          745586.7       90       127
## s2          712661.2       90        96
## Data attributes:
##       sim1             sim2             sim3             sim4       
##  Min.   : 7.827   Min.   : 7.458   Min.   : 7.789   Min.   : 7.346  
##  1st Qu.:12.930   1st Qu.:12.922   1st Qu.:12.890   1st Qu.:12.772  
##  Median :15.256   Median :15.129   Median :15.094   Median :15.176  
##  Mean   :15.844   Mean   :15.892   Mean   :15.670   Mean   :15.767  
##  3rd Qu.:18.217   3rd Qu.:18.520   3rd Qu.:17.967   3rd Qu.:18.220  
##  Max.   :32.247   Max.   :31.002   Max.   :29.252   Max.   :31.967

Symulacje geostatystyczne pozwalają również na przedstawianie niepewności interpolacji. W tym wypadku należy wykonać znacznie więcej powtórzeń (np. nsim = 100).

## drawing 100 GLS realisations of beta...
## [using conditional Gaussian simulation]

Uzyskane wyniki należy przeliczyć do oryginalnej jednostki, a następnie wyliczyć odchylenie standardowe ich wartości. Można to zrobić korzystając dwukrotnie z funkcji apply(). Przywrócenie oryginalnej jednostki odbywa się poprzez argumenty MARGIN = 2 oraz FUN = exp. Oznacza on, że funkcja exp() będzie wykonana na każdej kolumnie. W efekcie tej operacji liczba kolumn nie ulega zmianie.

Wyliczenie odchylenia standardowego odbywa się poprzez argumenty MARGIN = 1 oraz FUN = sd. W ten sposób funkcja sd() będzie wykonana na każdym wierszy. W efekcie tej operacji otrzymuje się tylko jedną kolumnę.

Finalnie otrzymujemy mapę odchylenia standardowego symulowanych wartości. Można na niej odczytać obszary o najpewniejszych (najmniej zmiennych) wartościach (niebieski kolor) oraz obszary o największej zmienności cechy (kolor żółty).

13.4 Sekwencyjna symulacja danych kodowanych

Symulacje geostatystyczne można również zastosować do danych binarnych. Dla potrzeb przykładu tworzona jest nowa zmienna temp_ind przyjmująca wartość TRUE dla pomiarów o wartościach temperatury niższych niż 12 stopni Celsjusza oraz FALSE dla pomiarów o wartościach temperatury równych lub wyższych niż 12 stopni Celsjusza.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.706  13.284  15.309  15.950  18.273  26.139
##    Mode   FALSE    TRUE 
## logical     212      38

W tej metodzie kolejne etapy przypominają przeprowadzenie krigingu danych kodowanych. Jedynie w funkcji krige() należy dodać argument mówiący o liczbie symulacji do przeprowadzenia (nsim).

##   model      psill    range
## 1   Nug 0.03299462    0.000
## 2   Sph 0.09866553 3006.406

## drawing 4 GLS realisations of beta...
## [using conditional indicator simulation]

Wynik symulacji danych kodowanych znacząco różni się od wyniku krigingu danych kodowanych. W przeciwieństwie do tej drugiej metody, w rezultacie symulacji nie otrzymujemy prawdopodobieństwa zajścia danej klasy, ale konkretne wartości 1 lub 0.

13.5 Zadania

  1. Stwórz nową siatkę dla obszaru o zasięgu x od 0 do 40000 i zasięgu y od 0 do 30000 oraz rozdzielczości 250. Wykonaj trzy symulacje bezwarunkowe
  2. Zbuduj po trzy symulacje bezwarunkowe w tej siatce korzystając z:
    • modelu sferycznego o semiwariancji cząstkowej 15 i zasięgu 7000
    • modelu nugetowego o wartości 1 razem z modelem sferycznym o semiwariancji cząstkowej 15 i zasięgu 7000
    • modelu sferycznego o semiwariancji cząstkowej 5 i zasięgu 7000
    • modelu sferycznego o semiwariancji cząstkowej 15 i zasięgu 1000

Porównaj graficznie uzyskane wyniki i opisz je. 3. Zbuduj po trzy symulacje bezwarunkowe w tej siatce korzystając z: - modelu potęgowego o semiwariancji cząstkowej 0.03 i zasięgu 0.3 - modelu potęgowego o semiwariancji cząstkowej 0.03 i zasięgu 0.6 - modelu potęgowego o semiwariancji cząstkowej 0.03 i zasięgu 0.9

Porównaj graficznie uzyskane wyniki i opisz je. 4. Stwórz optymalny model semiwariogramu zmiennej temp z obiektu punkty. Następnie korzystając z wybranej metody krigingowej, poznanej w poprzednich rozdziałach, wykonaj estymację zmiennej temp. 5. Wykonaj cztery symulacje warunkowe używając optymalnego modelu semiwariancji stworzonego w poprzednim zadaniu. Porównaj uzyskaną estymację geostatystyczną z symulacjami geostatystycznymi. Jakie można zaobserwować podobieństwa a jakie różnice? 6. Zbuduj optymalny model semiwariogramu empirycznego określający prawdopodobieństwo wystąpienia wartości ndvi poniżej 0.1 dla zbioru punkty. Wykonaj estymację metodą krigingu danych kodowanych. Następnie używając tego samego modelu, wykonaj cztery symulacje warunkowe (symulacje danych kodowanych). Jakie wartości progowe prawdobodobieństwa przypominają uzyskane symulacje?