LANL Earthquake Prediction – konkurs Kaggle

Od jakiegoś czasu już przymierzam się do napisania serii postów dotyczących konkursów uczenia maszynowego Kaggle. Czekałem na jakiś fajny konkurs, który będzie wdzięcznym tematem takich postów. Konkursy Kaggle mają to do siebie, że dotyczą na ogół bardzo wąskiego tematu i problemów z nim związanych. Są tam prawdziwe zbiory danych i prawdziwe problemy biznesowe. Ale najlepsze wyniki w takich konkursach uzyskiwane są przy pomocy modeli, które zdecydowanie nie nadają się jako samouczki, czy też nawet zaawansowane projekty end to end. Są po prostu bardzo złożone i niekoniecznie przejrzyste. Są również rozbudowywane metodą prób i błędów, więc wyglądają czasem jak potworki doktora Frankensteina. Czekałem więc na taki konkurs.

Ostatnio przeszło mi jednak przez myśl, że mogę się po prostu nie doczekać. Postanowiłem więc wziąć na tapet jeden z aktualnie otwartych konkursów. Nie wiem oczywiście, czy uda mi się cokolwiek sensownego w nim ugrać. Sądzę jednak, że pokazanie, że proces uczenia maszynowego to w większości czasu dość mizerne wyniki, może kogoś podbudować w jego dążeniach do rozwiązywania problemów.

Okej, koniec przynudzania, przyjrzyjmy się konkretom.

Parametry

• Nazwa: LANL Earthquake Prediction
• Opis: Can you predict upcoming laboratory earthquakes?
• Pula nagród: $50,000
• Sponsor: Los Alamos National Laboratory
• Zakończenie: 03.06.2019
• Metryka: mean absolute error
• Dane: 2GB (skompresowane)
• Typ danych: szeregi czasowe

Cel

Celem tego konkursu jest opracowanie metody na określenie, kiedy nastąpi trzęsienie ziemi, bazując na danych akustycznych z urządzenia pomiarowego. Dane te zostały przygotowane w laboratoriach LANL i współpracujących z nimi. Mamy więc szereg czasowy z informacją ile zostało do trzęsienia ziemi. Mamy też dane bez tej informacji – musimy ją wyznaczyć.

Mięso

Mięsem tego konkursu jest plik train.csv który waży 8,9G i wygląda tak:

acoustic_data,time_to_failure
12,1.4690999832
6,1.4690999821
8,1.469099981
5,1.4690999799
8,1.4690999788
8,1.4690999777
9,1.4690999766
7,1.4690999755
-5,1.4690999744

W pierwszej kolumnie mamy aktywność akustyczną (int16), a w drugiej dokładny czas do nastąpienia „trzęsienia ziemi” (float64). W pliku tym mamy 629_145_480 pomiarów. Są tam odczyty akustyczne poprzedzające wiele różnych trzęsień ziemi. I może mieć to chyba sens, bo z tego, co się orientuję, to nigdy nie ma po prostu jednego wstrząsu – jest to seria mniejszych i większych wstrząsów.

Dualnym dla pliku train.csv powinien być plik test.csv. W tym przypadku jest jednak inaczej. Dostajemy 2624 pliki o nazwach seg_*.csv, które wyglądają tak:

acoustic_data
5
2
9
12
6
5
2
0
5

Mamy więc tutaj tylko dane akustyczne – czego się spodziewaliśmy. Plik taki składa się z 150_000 pomiarów. Naszym zadaniem jest wyznaczyć czas do trzęsienia, dla każdego takiego pliku. Początek przykładowego pliku z rozwiązaniem wygląda tak:

seg_id,time_to_failure
seg_00030f,0
seg_0012b5,0
seg_00184e,0
seg_003339,0
seg_0042cc,0
seg_004314,0
seg_004cd2,0
seg_004ee5,0
seg_004f1f,0

Widzimy więc, że sama koncepcja przygotowania rozwiązania w tym konkursie będzie prosta.

Metryka

Metryką użytą w tym konkursie jest mean absolute error, czyli średni błąd bezwzględny. Jak się go oblicza? Załóżmy, że mamy trzy pliki do zbadania. Wyszło nam 3, 2 i 5 sekund. W rzeczywistości było to 3.2, 2 i 4. Mamy więc |3-3.2|, |2-2| i |5-4|. Średnia z tych błędów to (0.2+0+1)/3 = 0.4. Celem jest uzyskanie jak najniższej wartości. Widzimy więc, że ten element też jest dość intuicyjny w tym konkursie.

Benchmark

Jak to zwykle bywa z konkursami Kaggle, organizatorzy przygotowali benchmarki, które obrazują przykładowe rozwiązania na tablicy wyników. W tym konkursie mamy aktualnie trzy benchmarki:

• sample_submission.csv – 4.017 – wynik taki uzyskamy, jeśli jako rozwiązanie zgłosimy przykładowy plik z rozwiązaniem. Plik ten zawiera same 0, nie jest więc zbyt sensowny. Da się z niego natomiast wyliczyć średni błąd bezwzględny i tyle on właśnie wyniesie.
• Basic Feature Benchmark – 1.881 – jest to prosty model zaproponowany przez pracownika Kaggle. Bazuje na czterech cechach (min, max, mean i std) oraz na funkcji modelującej NuSVR z domyślnymi hiperparametrami. Kod tego rozwiązania omówię dalej w artykule.
• zixing AI – 1.481 – niestety nie udało mi się namierzyć kodu z tego rozwiązania.

Mamy więc tutaj niejako trzy progi wyznaczone przez konkretne benchmarki. Jak zaraz przekonasz się, środkowy benchmark jest całkiem prosty do samodzielnego osiągnięcia. Myślę więc, że śmiało można się zabierać za zaatakowanie poziomu zixing AI.

Basic Feature Benchmark

Przyjrzyjmy się więc środkowemu rozwiązaniu zaproponowanemu jako benchmark. Link do tego rozwiązania znajduje się tutaj. Zobaczmy komórka po komórce, co tutaj się dzieje:

# This Python 3 environment comes with many helpful analytics libraries installed
# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python
# For example, here's several helpful packages to load in 

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

# Input data files are available in the "../input/" directory.
# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory

import os
print(os.listdir("../input"))

# Any results you write to the current directory are saved as output.

['test', 'train.csv', 'sample_submission.csv']

Jest to standardowa pierwsza komórka, która jest w każdym nowym notebooku Kaggle. Mamy tutaj link do obrazu dockera, który jest używany. Mamy też podstawową informację o tym, jak używać notebooków, podstawowe importy oraz kod listujący zawartość katalogu input. Architekci tego systemu wymyślili sobie, że notebooki będą uruchamiane z katalogu workspace, dane będą trzymane w katalogu input, a wyniki w katalogu output. Jest to całkiem sensowne rozwiązanie, tylko trzeba się do niego przyzwyczaić ;-).

import matplotlib.pyplot as plt
from tqdm import tqdm
from sklearn.preprocessing import StandardScaler
from sklearn.svm import NuSVR
from sklearn.metrics import mean_absolute_error

Tutaj mamy dalsze importy modułów Pythonowych. W omawianym notebooku mamy więc:

• import numpy as np – podstawowy moduł numeryczny. Przydaje się tutaj np. do odpowiedniej definicji typów danych.
• import pandas as pd – moduł odpowiedzialny za ramki danych w Pythonie.
• import os – z tego, co widzę to jest to moduł użyty tutaj tylko do wypisywania listy plików.
• import matplotlib.pyplot as plt – interface do tworzenia wykresów.
• from tqdm import tqdm – graficzna nakładka pokazująca postęp pętli for.
• from sklearn.preprocessing import StandardScaler – funkcja skalująca od -1 do 1.
• from sklearn.svm import NuSVR – funkcja obsługująca Nu Support Vector Regression.
• from sklearn.metrics import mean_absolute_error – funkcja wyliczająca naszą docelową metrykę.

Jak widać, dużo tego nie ma i nie są to też egzotyczne moduły.

train = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})

Tutaj mamy standardowe wczytywanie do ramki danych z pliku csv. Warto zwrócić uwagę na to, że jeżeli znamy typy danych, na których będziemy operować, to możemy je tutaj od razu podać. W ten sposób Pandas nie będzie zgadywać ile pamięci zarezerwować dla przyszłej ramki, tylko od razu skorzysta z naszych wskazań. Jeśli wiemy, co robimy, to minimalizujemy ryzyko nadmiernego zużycia pamięci RAM.

train.head()

	acoustic_data	time_to_failure
0	12	1.4691
1	6	1.4691
2	8	1.4691
3	5	1.4691
4	8	1.4691

Tutaj podglądamy sobie początek naszej ramki danych.

# pandas doesn't show us all the decimals
pd.options.display.precision = 15

Okazuje się, że uzyskane wyjście nie pokazuje liczb z taką dokładnością, jaką mamy w ramce danych. Przekazujemy więc do Pandas odpowiednio zmieniony parametr odpowiedzialny za wyświetlanie wartości „po przecinku”.

# much better!
train.head()

	acoustic_data	time_to_failure
0	12	1.4690999832
1	6	1.4690999821
2	8	1.4690999810
3	5	1.4690999799
4	8	1.4690999788

Teraz już wyświetlone liczby zgadzają się z zawartością pliku. Wyznaczmy więc powyżej wspomniane cechy:

# Create a training file with simple derived features

rows = 150_000
segments = int(np.floor(train.shape[0] / rows))

X_train = pd.DataFrame(index=range(segments), dtype=np.float64,
                       columns=['ave', 'std', 'max', 'min'])
y_train = pd.DataFrame(index=range(segments), dtype=np.float64,
                       columns=['time_to_failure'])

for segment in tqdm(range(segments)):
    seg = train.iloc[segment*rows:segment*rows+rows]
    x = seg['acoustic_data'].values
    y = seg['time_to_failure'].values[-1]
    
    y_train.loc[segment, 'time_to_failure'] = y
    
    X_train.loc[segment, 'ave'] = x.mean()
    X_train.loc[segment, 'std'] = x.std()
    X_train.loc[segment, 'max'] = x.max()
    X_train.loc[segment, 'min'] = x.min()

100%|██████████| 4194/4194 [00:08<00:00, 469.57it/s]

Widzimy, że kluczowym elementem tej komórki jest pętla for idąca po segmentach. O co tutaj chodzi? Otóż jeśli wrócimy pamięcią do opisu plików testowych, to przypomnimy sobie, że pliki te nazywają się seg_ i jakiś id. Każdy z tych plików ma 150_000 odczytów. Chodzi więc o to, że jeśli będziemy dokonywać predykcji, to zawsze będziemy to robić dla jednego pliku zawierającego 150_000 wartości. A nasz plik treningowy ma tych wartości miliony. Powinniśmy więc w momencie trenowania doprowadzić do podobnej sytuacji – funkcja modelująca powinna widzieć wiele „plików zawierających 150_00 wartości. W zmiennej segments trzymamy więc liczbę możliwych do uzyskania rozłącznych segmentów (nieco ponad 4 tysiące w tym przypadku). Przechodzimy więc przez duży plik train, wyciągamy z niego segment i wyliczamy dla niego pożądane cechy – mean, std, max i min. W ten sposób porzucamy niejako oryginalne pełne dane i przechodzimy na segment i cztery cechy z nim związane. Pozostaje nam jeszcze zmienna zależna. Tutaj bierzemy po prostu ostatnią wartość z segmentu.

X_train.head()

	ave	std	max	min
0	4.884113333333334	5.101089126891323	104.0	-98.0
1	4.725766666666667	6.588801819164257	181.0	-154.0
2	4.906393333333333	6.967373808828945	140.0	-106.0
3	4.902240000000000	6.922282112791032	197.0	-199.0
4	4.908720000000000	7.301085852684289	145.0	-126.0

Widzimy tutaj, że mamy teraz nieco inną sytuację niż początkowo. Wcześniej mieliśmy szereg czasowy z przyporządkowaną do niego liczbą, a teraz mamy obserwacje z czterema cechami. Dość mocno uprościliśmy sytuację poprzez usuwanie informacji.

scaler = StandardScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)

Teraz autor zaproponował wyskalowanie uzyskanych cech do przedziału -1 i 1. I to by było na tyle, jeśli chodzi o przygotowanie cech na bazie danych treningowych. Faktyczny trening wygląda standardowo:

svm = NuSVR()
svm.fit(X_train_scaled, y_train.values.flatten())
y_pred = svm.predict(X_train_scaled)

Mamy tutaj konstruowanie modelu z domyślnymi hiperparametrami, jego trening i dokonywanie predykcji. Możemy wziąć wszystkie te przewidziane wartości i porównać je z prawdziwymi wartościami:

plt.figure(figsize=(6, 6))
plt.scatter(y_train.values.flatten(), y_pred)
plt.xlim(0, 20)
plt.ylim(0, 20)
plt.xlabel('actual', fontsize=12)
plt.ylabel('predicted', fontsize=12)
plt.plot([(0, 0), (20, 20)], [(0, 0), (20, 20)])
plt.show()

Oprócz tego możemy też wyliczyć sobie score z wykorzystaniem metryki konkursowej:

score = mean_absolute_error(y_train.values.flatten(), y_pred)
print(f'Score: {score:0.3f}')

Score: 2.314

Pozostaje nam już tylko dokonanie finalnych predykcji dla plików testowych dostarczonych przez Kaggle. Najłatwiej będzie zacząć od wczytania pliku sample_submission.csv w celu systematycznego uzupełnienia go naszymi wynikami.

submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')

Musimy sobie jeszcze utworzyć ramkę danych z danymi, które wrzucimy w wytrenowany model. Musi mieć ona tyle samo kolumn co ramka użyta w treningu, możemy więc śmiało się o nią oprzeć. Indeksem w tej ramce będzie ten sam indeks co w powyższym pliku sample_submission.csv.

X_test = pd.DataFrame(columns=X_train.columns, dtype=np.float64, index=submission.index)

Teraz czeka nas najtrudniejsza część przygotowania pliku wynikowego. Musimy przejść po każdym pliku testowym i wygenerować dla niego odpowiednie cechy. Najlepiej zrobić to w pętli for idącej po indeksie ramki danych, którą będziemy uzupełniać. Indeks ten to również nazwy plików, dla każdego wiersza wczytujemy więc plik csv i na jego podstawie generujemy cztery cechy, które umieszczamy w odpowiednich kolumnach. Wygląda to tak:

for seg_id in X_test.index:
    seg = pd.read_csv('../input/test/' + seg_id + '.csv')
    
    x = seg['acoustic_data'].values
    
    X_test.loc[seg_id, 'ave'] = x.mean()
    X_test.loc[seg_id, 'std'] = x.std()
    X_test.loc[seg_id, 'max'] = x.max()
    X_test.loc[seg_id, 'min'] = x.min()

Jak już mamy wyznaczone cechy dla danych testowych to możemy przejść do dokonania predykcji i zapisania wyników w pliku na dysku. Nie możemy jednak zapomnieć o zastosowaniu tej samej funkcji skalującej co wcześniej (scaler.transform(X_test)). Bez tego przekażemy dane z zupełnie innymi rzędami wielkości i nie będzie mało to za wiele sensu.

X_test_scaled = scaler.transform(X_test)
submission['time_to_failure'] = svm.predict(X_test_scaled)
submission.to_csv('submission.csv')

Notebook z dokładnie taką zawartością da wynikowy plik submission.csv, który po wysłaniu do konkursu da wynik 1.881. Możesz też od razu uruchomić ten benchmark na platformie Kaggle, klikając „Fork” a później „Commit” na stronie tego notebooka tutaj.

Potencjalne ulepszenia powyższego benchmarku

Powyższy benchmark pozostawia wiele pola do popisu. Możemy tutaj podjąć co najmniej trzy działania, żeby znaleźć lepszy wynik:

• optymalizacja hiperparametrów – jak widzieliśmy, nie było tutaj żadnego wyboru hiperparametrów. Bierzemy funkcję modelującą z domyślnymi wartościami hiperparametrów i tyle. Jeśli postanowimy dobrać lepsze wartości hiperparametrów może uda się podkręcić trochę ten wynik.
• inna funkcja modelująca – skorzystaliśmy tutaj z funkcji NuSVR. A jak wiemy, sam Scikit-Learn oferuje ich całą masę. Nie mówiąc już o innych modułach albo własnym kodzie. Może więc warto sprawdzić inne funkcje?
• cechy – stworzenie nowych cech. Cechy, które tutaj mamy są dość ubogie. Dodanie nowych cech albo całkiem zmienienie sposobu ich tworzenia może dać lepsze rezultaty.

Oczywiście nie trzeba wcale bazować na tym benchmarku. Można stworzyć rozwiązanie w dowolny sposób. A z racji tego, że Kaggle w swoim środowisku udostępnianym na potrzeby tego i innych konkursów oferuje nawet dostęp do GPU (info tutaj – sekcja Technical Specifications), to można się pokusić nawet o deep learning. Mnie, na razie, nie udało się uzyskać znacząco lepszego wyniku, ale jeśli taki uzyskam, to nie omieszkam podzielić się nim na tym blogu.

Jeśli interesuje Cię jakiś temat – nie musi być związany z tym artykułem – to zostaw mi sygnał tutaj. Dzięki!