Przewidywanie pogody

W poprzednim odcinku zebraliśmy dane o pogodzie w Polsce od 1951 roku. Dzisiaj przygotujemy prognozę pogody. Ale przede wszystkim nauczymy się weryfikować modele.

Mając zgromadzone dane o pogodzie z ponad 60 lat (1951-2017) możemy pokusić się o prognozowanie temperatury w przyszłości. Wykorzystamy trzy (właściwie dwa) modele dostępne od ręki w popularnych bibliotekach R oraz zaprzęgniemy do prognozowania sieć neuronową.

Zaczniemy od potrzebnych bibliotek:

Na początek pobieramy dane o średniej dziennej temperaturze wszystkich stacji w latach 2000-2016. Te dane potraktujemy jako dane treningowe. Dane z 2017 roku potraktujemy jako testowe.

Na potrzeby zbudowania modeli ARIMA i ETS potrzebujemy danych w postaci szeregu czasowego. Budujemy taki z danych treningowych:

Możemy szereg czasowy w bardzo prosty sposób obejrzeć i z grubsza przeanalizować składowe, między innymi sezonowość i trend:

Pięknie widać sezonowość temperatury. Widać też, że w okolicy 2006 roku było zdecydowanie cieplej (składowa trendu). Można to porównać z wykresami z poprzedniej części; na tym wykresie:

oraz na tym:

Spróbujmy przewidzieć temperaturę na przyszłość, korzystając z danych historycznych. Nie jest to może najmądrzejsze podejście do prognozowania pogody, ale skoro temperaturę charakteryzuje sezonowość, to dlaczego nie?

Model ARIMA (auto)

Pierwszy będzie model typu ARIMA z automatycznie dobranymi współczynnikami. O modelu tym pisałem jakiś rok temu, przy okazji prognozowania stopy bezrobocia, który to model sprawdził się bardzo dobrze (porównajcie ówczesne prognozy z danymi o stopie bezrobocia w 2017 roku).

Model nie dał żadnych sensownych wyników – temperatura na początku roku urośnie, a potem będzie już stała. W naszym klimacie to bzdura (chociaż dzisiaj klikałem po nowej wersji serwisu z moją ulubioną prognozą – meteo.pl – i na przykład we Włoszech już teraz w nocy jest tak samo ciepło jak w dzień).

Model ARIMA

Należałoby dostosować parametry modelu, przede wszystkim parametr mówiący o długości okresu (sezonu) – w naszym przypadku jest to rok (365.25 dni).

Jest zdecydowanie lepiej niż w przypadku automatycznego modelu ARIMA, widać już górkę w lecie i spadek w zimie. Linia niebieska to prognoza (i tego użyjemy do dalszych porównań), a szarości to przestrzeń w której teoretycznie temperatura może się poruszać z określonym stopniem ufności.

Model ETS

W prognozowaniu bezrobocia świetnie sprawdził się model ETS, ale dla tak dużej częstotliwości szeregu czasowego przy wywołaniu funkcji ets() dostajemy komunikat, że to jednak przesada (z tą częstotliwością) i należy użyć stlf(). No to proszę bardzo:

Na pierwszy rzut oka też jest dobrze. Co ciekawe – prognozy przeliczyły się zdecydowanie szybciej niż paramteryzowana Arima.

Przygotujmy sobie tabelkę z wartościami prognozowanymi przez kolejne modele – przyda się za chwilę do porównania który model daje najlepsze predykcje.

Sieć LSTM

Ostatnim modelem będzie model zbudowany w oparciu o sieć neuronową typu LSTM. Zadziałamy podobnie jak w przypadku przewidywania tekstu: podzielimy historię na serię okien (tutaj o długości 30 dni) i wytrenujemy sieć tak, aby nauczyła się odpowiadać na wektor 30-elementowy jedną liczbą (temperaturą z kolejnego dnia).

Pierwszy krok to przygotowanie danych w odpowiedniej dla TensorFlow (w tym przypadku opakowany w bibliotekę Keras) formie:

Teraz czas na stosowny model. Próbowałem kilku architektur sieci, ta daje dobre wyniki przy niewielkim czasie obliczeń.

Podsumowanie modelu (architektura sieci) wygląda następująco:

Mamy prawie 400 tysięcy parametrów do określenia. Zatem do dzieła, licz komputerku, licz!

W przypadku sieci zawsze warto sprawdzić jak wyglądają krzywe błędów:

Wygląda to dobrze, krzywe spadają z liczą epochów, a w dodatku są ze sobą zbieżne – nie ma więc mowy o przeuczeniu sieci. Możecie popróbować zmienić liczbę epochów oraz węzłów sieci.

Policzmy więc temperaturę dla 2017 roku. Przygotujmy dane w odpowiedniej formie:

i sprawdźmy wynik:

Teraz do zestawienia wszystkich prognoz dodajemy te z modelu LSTM. Ważna sprawa – okno 30-dniowe sprawia, że nie mamy pełnego roku, zatem aby nie było problemów część danych z poprzednich modeli sobie po prostu usuniemy. To szkoła jest, a nie poważny instytut meteo, więc możemy :)

Porównanie wyników

Zobaczmy jak według poszczególnych modeli miała wyglądać temperatura w 2017 roku. Oczywiście w zestawieniu z tą prawdziwą – inaczej test nie miałby większego sensu.

Czerwona linia to rzeczywiste średnie dzienne temperatury z 2017, punkty – temperatury według kolejnych modeli.

Automatyczna Arima dała ciała, co wiemy już od początku. Cała reszta leży mniej więcej w okolicach rzeczywistości. Ciekawe są czerwone punkty – odpowiadające za model LSTM. Wyglądają trochę jak średnia ruchoma.

Ale miarą poprawności nie jest mniej więcej a rzeczywiste różnice od prawdziwych wartości. Matematycznie zwie się to residua. Policzmy residua, a dodatkowo na jednym wykresie pokażmy jak układają się wartości obliczone przez model w porównaniu z wartościami rzeczywistymi.

Na powyższym wykresie interesuje nas tak na prawdę pierwszy wiersz (albo pierwsza kolumna). To są wykresy typu Q-Q plot – w idealnym modelu punkty ułożyłyby się na linii prostej (takiej na ukos, z równaniem y=x). Widać, że najmniejszy rozrzut wokół prostej jest dla pary org – lstm, co może prowadzić do wniosku, że model LSTM dał najmniejsze błędy.

RMSE

Czy tak jest również można policzyć – używając na przykład miary RMSE (root mean square error – błąd średnio kwadratowy), czyli pierwiastka z sumy kwadratów odchyłek (residuów). RMSE dla kolejnych modeli to:

  • ARIMA (auto): 133.4418
  • ARIMA (parametryzowana): 85.377
  • STLF: 73.9815
  • LSTM: 50.6776

Mamy potwierdzenie, że LSTM dał najmniejsze odchyłki.

Można to też sprawdzić poprzez współczynniki korelacji danych obliczonych z rzeczywistymi – w tabelce:

org arima arima_params stlf lstm
org 1.000 0.716 0.822 0.878 0.948
arima 0.716 1.000 0.552 0.639 0.693
arima_params 0.822 0.552 1.000 0.919 0.850
stlf 0.878 0.639 0.919 1.000 0.917
lstm 0.948 0.693 0.850 0.917 1.000

albo na wykresie:

Znowu – interesuje na pierwszy wiersz lub pierwsza kolumna. Widzimy, że współczynnik korelacji pomiędzy org a lstm jest największy. Bardzo fajnie widać to na miniaturkach (obrazują to co pairs() wyżej).

Różnice dzień po dniu

Na koniec sprawdźmy jak bardzo dla każdego dnia modele rozbiegały się z rzeczywistością.

Warto zwrócić uwagę na kierunek słupków na poszczególnych wykresach w konkretnych dniach (albo okresach) – z reguły modele mylą się w tą samą stronę – czyli albo wszystkie prognozują temperaturę za wysoką w porównaniu z rzeczywistością, albo za niską. To daje jakieś spojrzenie na anomalia temperatury (rzeczywistej) w danym dniu na przestrzeni wielu lat.

Rozkład residuów

Sprawdźmy czy różnice rozkładają się równomiernie na plus i na minus.

Im bardziej smukła krzywa i jej wierchołek bliżej zera tym lepiej. Bardzo dobrze wygląda tutaj rozkład residuów dla modelu lstm – po raz kolejny wygrywa. Zielona arima_res już na tym wykresie wygląda nieciekawie i możemy sądzić, że jest nieprzydatna. Zweryfikujmy to przy pomocy testu t-Studenta.

Zazwyczaj dwie średnie z różnych od siebie grup będą się różnić. Test t-Studenta powie nam czy te różnice są istotne statystycznie. O co chodzi? Dla tych, którzy niewiele pamiętają ze statystyki prostu przykład – czy kobiety i mężczyźni palą tyle samo (średnio)? Stawiamy dwie przeciwne hipotezy – H0 oraz H1.

Hipoteza zerowa takiego testu będzie brzmiała:

  • H0: Średnia liczba wypalanych papierosów w grupie mężczyzn jest taka sama jak średnia liczba wypalanych papierosów w grupie kobiet.

Hipoteza alternatywna z kolei jest przeciwna:

  • H1: Kobiety będą różnić się od mężczyzn pod względem liczby wypalanych papierosów w ciągu miesiąca (licząc średnią).

Jeśli wynik testu t-Studenta będzie istotny na poziomie p < 0.05 możemy odrzucić hipotezę zerową na rzecz hipotezy alternatywnej (H0 jest nieprawdziwe, a z tego wynika, że H1 jest prawdziwe). Czyli mężczyźni różnią się istotnie od kobiet jeśli chodzi o liczbę wypalanych papierosów. Patrząc na średnie w obu grupach dowiemy się kto pali więcej.

Wróćmy do naszych danych o pogodzie. Porównamy średnie temperatury uzyskane różnymi modelami ze średnią temperaturą zarejestrowaną rzeczywiście. Tutaj porównujemy dane (średnie) z całego roku, ale warto rozbić te porównania na poszczególne pory roku albo miesiące. Później na takich grupach policzyć na przykład ile razy hipoteza H0 została odrzucona dla danego modelu. Tam, gdzie ta wartość (liczba grup z odrzuconą hipotezą H0) jest większa mamy gorszy model.

Spróbujmy najpierw rozpatrzyć średnią z temperatur dla całego roku, bez dzielenia na grupy:

ARIMA (auto):

Wartość p jest mniejsza od 0.05, a więc odrzucamy H0 (średnia roczna temperatura jest różna dla danych z modelu i danych rzeczywistych).

  • ARIMA (parametryzowana): p = 0.9791, a więc nie odrzucamy H0
  • STLF: p = 0.069, a więc nie odrzucamy H0
  • LSTM: p = 0.2315, a więc nie odrzucamy H0

Model LSTM wydaje się być jak na razie najlepszy, gdyż:

  • nie możemy odrzucić hipotezy zerowej H0, a więc średnia roczna według modelu nie różni się statystycznie od średniej zaobserwowanej
  • model ten ma największy współczynnik korelacji z danymi rzeczywistymi (co łączy się bezpośrednio z punktem wyżej – być może nawet jest na to jakiś dowód matematyczny, albo jedno wynika z drugiego?)
  • model ten ma najniższą wartość RMSE

Na (kolejny już) pierwszy rzut oka jest zwycięzca. Spróbujmy jeszcze przeprowadzić test t-Studenta na grupach dla poszczególnych miesięcy.

Mamy 6 miesięcy (z 11) dla których hipoteza o równości średnich (miesięcznych) z modelu LSTM i rzeczywistych jest prawdziwa. Tutaj wynik jest różny i może być losowy. Warto ustawić na początek generator liczb losowych na stałą wartość (przez set.seed())

Dla pozostałych modeli H0 jest prawdziwe w przypadku

  • ARIMA (auto): 3 miesięcy
  • ARIMA (parametryzowana): 4 miesięcy
  • STLF: 5 miesięcy

Znowu LSTM wydaje się być najlepszym.

Pięknie, pięknie ale jest tutaj jeden poważny błąd. Jaki?

Ano taki, ze przewidujemy wartości w LSTM korzystając z wartości prawdziwych. Bierzemy wektor (jedno okno) wartości przewdziwych, obliczamy jedną wartość przewidywaną. Później bierzemy kolejne okno i znowu przewidujemy jedną wartość. Zapamiętujemy wartości przewidziane na podstawie wartości prawdziwych. Robimy więc trochę taką prognozę krótkoterminową i tę prognozę porównujemy z rzeczywistością.

W przypadku modeli zbudowanych na szzeregach czasowych (ARIMA i STFL) przewidujemy temperaturę na cały 2017 rok od razu. Czy zatem można porównywać w ten sposób modele?

Może w modelu LSTM trzeba przewidywać kolejne dni na podstawie już przewidzianych? To brzmi sensowniej. Trochę gimnastyki i z tym sobie poradzimy:

Uuuu – trochę bullshit! Nie widać lata (powinno być w górę w okolicach 100-200, później znowu w dół). W przewidywaniu tekstu mięliśmy podobny problem, pomogła zmiana architektury sieci. Tutaj problem może wynikać z tego, że za początkowe (lstm_predict ustalane na początku) dane bierzemy temperatury z grudnia 2016 (gdzie raczej robiło się coraz chłodniej – stąd coraz niższe przewidywania). Być może zaczynając od czerwca wynik byłby odwrotny. Pomóc też może wcześniejsze (przed trenowaniem) skalowanie danych. Dla zainteresowanych polecam post Using LSTMs to forecast time-series – tam znajdziecie podejście do tematu z przykładowym kodem w Pythonie. A może po prostu nie można prognozować tak daleko do przodu?

Na zakończenie podsumujmy, czego nauczyliśmy się dzisiaj:

  • Do prognozowania szeregów czasowych przydatny jest pakiet forecast
  • Po przygotowaniu predykcji warto przeanalizować wyniki na próbce testowej – szczególnie residua, ich rozkład. Przydatne mogą być współczynniki korelacji i test t-Studenta. I to jest najważniejsze w dzisiejszym poście.
  • Sieć neuronowa typu LSTM też może być przydatna, szczególnie jeśli przewidujemy krótkie okresy (a nie cały rok)

Uwagi? Komentarze? Zapraszam do dyskusji poniżej oraz na fanpage Dane i analizy.

1 komentarz do wpisu “Przewidywanie pogody

Dodaj komentarz