{"cells":[{"metadata":{},"cell_type":"markdown","source":"Projekt rozpoczeliśmy od podstawowych czynności takich jak zdefiniowanie ścieżek do wczytywania plików, podstawowe importy takich bibliotek jak numpy oraz pandas. Początkowo wczytaliśmy tylko 1/2 danych ze zbioru treningowego z powodu sporych rozmiarów jak i problemów wydajnościowych. Rozmiar zbioru treningowego przekracza 10GB, co uniemożliwiało normalną pracę na pliku. Połowa zbioru jest wystarczająca do zaprezentowania podstawowych miar tendecji centralnych, z których później byliśmy w stanie wyciągnąć odpowiednie wnioski. Dodatkowo wczytaliśmy dane ze zbioru testowego oraz sprawdziliśmy ile próbek zawierał każdy segment danych testowych. Każdy z plików zawierał 150000 próbek. Warto również wspomnieć, że dane treningowe oraz testowe, które zostały udostępnione na potrzeby konkursu bardzo od siebie odbiegają, gdyż pomiary dla danych testowych nie są bezpośrednią kontynuacją pomiaru dla danych treningowych. \n\nWyniki eksploracji przedstawiliśmy w kolejnych etapach prac na wykresach oraz tabelach"},{"metadata":{"trusted":false},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np\nimport pandas as pd \n\nimport os\n\nDATA_DIR = \"../input\"\nTEST_DIR = r'../input/test'\n\ntrain_df = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'), nrows=300000000, dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\n\nprint(train_df.info())\n\n\nld = os.listdir(TEST_DIR)\nsizes = np.zeros(len(ld))\n\nfor i, f in enumerate(ld):\n    df = pd.read_csv(os.path.join(TEST_DIR, f))\n    sizes[i] = df.shape[0]\n\nprint(np.mean(sizes))\nprint(np.min(sizes))\nprint(np.max(sizes))\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Poniżej zebraliśmy wszystkie potrzebne i wykorzystywane biblioteki oraz frameworki, które ułatwiały nam pracę w obliczeniach, rysowaniu wykresów, klasyfikacji, regresji, eksploracji i analizie danych, predykcji.\n\nBiblioteki Seaborn oraz Matplotlib wykorzystaliśmy do graficznej wizualizacji wykresów.\nBiblioteka Catboost znalazła zastosowanie w naszym projekcie w zadaniach związanych z regresją, uczeniem maszynowym.\nNajszerszym zakresem zastosowań okazały się bilioteki sklearn, numpy i seaborn do uczenia maszynowego, wizualizacji i operacji matematycznych, które zapewniła świetne narzędzia do eksploracji i analizy danych."},{"metadata":{"trusted":false},"cell_type":"code","source":"import seaborn as sns\nimport matplotlib.pyplot as plt\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV, KFold\nfrom sklearn.svm import SVR, NuSVR\nfrom sklearn.kernel_ridge import KernelRidge\nimport pandas as pd\nimport numpy as np\nimport os\nimport gc\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nDATA_DIR = \"../input\"\nTEST_DIR = r'../input/test'\n\nld = os.listdir(TEST_DIR)\nsizes = np.zeros(len(ld))\n\nfrom scipy.signal import hilbert\nfrom scipy.signal import hann\nfrom scipy.signal import convolve\nfrom scipy.stats import pearsonr\nfrom scipy import stats\nfrom sklearn.kernel_ridge import KernelRidge\n\nimport lightgbm as lgb\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_absolute_error\n\nfrom tsfresh.feature_extraction import feature_calculators\n\n%matplotlib inline\nsns.set_style('darkgrid')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Początkowo rysowaliśmy wykresy z połowy zbioru treningowego, a w tym konkretnym przypadku poniżej braliśmy co 50 daną z zbioru treningowego. Na czerwono zaznaczone zostały dane akustyczne(acoustic_data), a na niebiesko czas(time_to_failure). Pomimo tego wciąż mogliśmy zauważyć pewne zależności. Na wykresie poniżej przedstawiliśmy zależność danych akustycznych(acoustic_data) od czasu trzęsienia(time_to_failure). Są to dane w czystej postaci - jeszcze niemodyfikowane na tym etapie. Jak widać przy każdym szczycie(peak) danych akustycznych widzimy momentalny zryw(peak) czasu od wartości niemalże zerowych do góry, który symbolizuje trzęsienie ziemi, a następnie liniową zmianę czasu biegnącą ponownie do zera. \n\nTak jak widać na poniższym obrazku można znaleźć 7 takich sytuacji, które odpowiadają 7 trzęsieniom w tej konkretnej wybranej części zbioru treningowego. Co najistotniejsze w tym wykresie to fakt, że szczyt wartości sygnału występuje przed prawie każdym wystąpieniem trzęsienia. Z tego faktu na pewno warto będzie zwrócić uwagę na wykresy przedstawiające wartości maksymalne oraz odchylenie standardowe.\n\nWarto wspomnieć, że w całym zbiorze treningowych samych trzęsień jest tylko 16, w stosunku do wszystkich pomiarów jest to wartośc poniżej 1% !"},{"metadata":{"trusted":false},"cell_type":"code","source":"acoustic_data_sample = train_df['acoustic_data'].values[::50]\ntime_to_failure_sample = train_df['time_to_failure'].values[::50]\n\nfig, ax1 = plt.subplots(figsize = (12,8))\nplt.title('Data from DF')\nplt.plot(acoustic_data_sample, color='r')\nax1.set_ylabel('Acousting Data', color='r')\nplt.legend(['acoustic data'], loc=(0.01, 0.95))\nax2 = ax1.twinx()\nplt.plot(time_to_failure_sample, color='b')\nax2.set_ylabel('Time to Failure', color='b')\nplt.legend(['time_to_failure'], loc=(0.01, 0.95))\nplt.grid(True)\n\ndel acoustic_data_sample, time_to_failure_sample\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Tak jak wspomnieliśmy wyżej na wykresie widoczne jest 7 trzęsień z wszystkich 16, które występują w pełnym zbiorze danych. Założyliśmy, że praca na danych, w który występuje prawie połowa trzęsień może przynieść podobne rezultaty, a będzie na pewno szybsza w porównaniu do pracy na pełynm zbiorze. Analizując tylko część zbioru możemy na pewno natrafić później na problemy związane z faktem, że po każdym skoku wartości danych akustycznych występuje trzęsienie. Rozważyliśmy dodanie cechy(feature), która zawierałaby wartość czasową od momentu nastąpienia szczytu(peak) wartości akustycznej, a faktycznego trzęsienia. Czasy dla wszystkich byłyby bardzo zblizone i analizując kolejne wykresy wynosiłyby średnio 0.31 ms."},{"metadata":{},"cell_type":"markdown","source":"Następnie wygenerowaliśmy losowe ziarno oraz tablice losowych wartości z przedziału odpowiadającemu połowie zbioru treningowego. Rysowaliśmy wykres z krokiem co 150000, ponieważ tyle wynosiła ilość próbek w każdym z pomiarów. Jeśli przyjrzeć się 9 poniższym wykresom możemy zauważyć, że wyniki danych akustycznych w większości przypadków z zakresu (-200;300). Jednak jak widać wystąpiły też pomiary z wynikami o wiele większymi np.: wynik w szczycie koło 1500 oraz około 4000."},{"metadata":{"trusted":false},"cell_type":"code","source":"np.random.seed(2018)\nrand_idxs = np.random.randint(0, 300000000-150000, size=9, dtype=np.int32)\nfig, axes = plt.subplots(3, 3, figsize=(18, 10))\n\nfor x in range(3):\n    for y in range(3):\n        ad = train_df['acoustic_data'].values[rand_idxs[x*3 + y]: rand_idxs[x*3 + y] + 150000]\n        ttf = train_df['time_to_failure'].values[rand_idxs[x*3 + y]: rand_idxs[x*3 + y] + 150000]\n\n        axes[x][y].plot(ad, color='blue')\n        axes[x][y].set_xticks([])\n\n        s = axes[x][y].twinx()\n        s.plot(ttf, color='red')\n        \nplt.tight_layout()\nplt.show()\ndel ad, ttf, rand_idxs\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"W tym miejscu wyznaczyliśmy sobie jaka jest liczba trzęsień ziemi w danych, na których pracujemy oraz liczba szczytów(peak) w tychże danych. Liczby sobie odpowiadały także wyznaczyliśmy indeksy tych szczytów."},{"metadata":{"trusted":false},"cell_type":"code","source":"d = {'vals': train_df[train_df['acoustic_data']>2000].index.values}\npeaks = pd.DataFrame(d)\npeaks['diff'] = peaks['vals'].diff(periods=-1)\nselected_peaks = peaks[abs(peaks['diff'])>30000]['vals'].values\n\n\ntrain_df['diff'] = train_df['time_to_failure'].diff()\nindexes_of_eartgquakes = train_df[train_df['diff']>1].index.values\n\nprint('Number of earthquakes in loaded data: ', len(indexes_of_eartgquakes))\nprint('Number of peaks in loaded data: ',len(selected_peaks))\n\nprint(selected_peaks)\n\ndel peaks, d\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Po wyznaczeniu interesujących nas indeksów, narysowaliśmy wykresy tych wartości w próbce 150000, aby dokładniej przyjrzeć się znalezionym szczytom. Wspominany już wielokrotnie szczyt sygnału tuż przed trzęsieniem powtarza się dosyć często. Jest to na pewno ważna i przydatna infmoracja, jeśli jednak weźmiemy pod uwagę, że takich sytuacji w całym zbiorze podzielonym na kawałki z próbkami liczącymi 150000 jest jedynie 16, może okazać się, że ciężko będzie ją wykorzystać w odpowiedni sposób na reszcie danych."},{"metadata":{"trusted":false},"cell_type":"code","source":"fig, axes = plt.subplots(4, 2, figsize=(18, 10))\nfig.delaxes(axes[3,1])\n\nfor x in range(7):\n        ad = train_df['acoustic_data'].values[selected_peaks[x]-75000: selected_peaks[x]+75000]\n        ttf = train_df['time_to_failure'].values[selected_peaks[x]-75000: selected_peaks[x]+75000]\n\n        axes[int(x/2)][x%2].plot(ad, color='blue')\n        axes[int(x/2)][x%2].set_xticks([])\n\n        s = axes[int(x/2)][x%2].twinx()\n        s.plot(ttf, color='red')\n        \nplt.tight_layout()\nplt.show()\ndel ad, ttf\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Poniżej wykresy przedstawiające dane treningowe w dziedzine czasu, na których widać momenty tręsień ziemi, z których niestety cięzko jest cokolwiek wywnioskować."},{"metadata":{"trusted":false},"cell_type":"code","source":"fig, axes = plt.subplots(4, 2, figsize=(18, 10))\nfig.delaxes(axes[3,1])\n\nfor x in range(7):\n        ad = train_df['acoustic_data'].values[indexes_of_eartgquakes[x]-140000: indexes_of_eartgquakes[x]+10000]\n        ttf = train_df['time_to_failure'].values[indexes_of_eartgquakes[x]-140000: indexes_of_eartgquakes[x]+10000]\n\n        axes[int(x/2)][x%2].plot(ad, color='blue')\n        axes[int(x/2)][x%2].set_xticks([])\n\n        s = axes[int(x/2)][x%2].twinx()\n        s.plot(ttf, color='red')\n        \nplt.tight_layout()\nplt.show()\ndel ad, ttf\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### W dalszej części analiz przystapimy do wyznaczania feature'ów do naszego przyszłego modelu"},{"metadata":{},"cell_type":"markdown","source":"Algorytm STA/LTA(Short Time Average over Long Time Average) zaprojektowany jest w celu ignorowania sygnałów energetycznych z pomiarów wibracji otoczenia. Wybiera on tylko te części sygnałów, których amplituda jest pseudo-stacjonarna. Algorytm zmierza do uniknięcia wybuchów energii. Poniżej znajduje się jego implementacja."},{"metadata":{"trusted":false},"cell_type":"code","source":"def classic_sta_lta(x, length_sta, length_lta):\n    \n    sta = np.cumsum(x ** 2)\n\n    # Zamiana na float\n    sta = np.require(sta, dtype=np.float)\n\n    # Kopia dla LTA\n    lta = sta.copy()\n\n    # Obliczanie STA i LTA\n    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n    sta /= length_sta\n    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n    lta /= length_lta\n\n    # Uzupełnienie zerami\n    sta[:length_lta - 1] = 0\n\n    # Aby nie dzielić przez 0 ustawiamy 0 na małe liczby typu float\n    dtiny = np.finfo(0.0).tiny\n    idx = lta < dtiny\n    lta[idx] = dtiny\n\n    return sta / lta","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def calc_change_rate(x):\n    change = (np.diff(x) / x[:-1]).values\n    change = change[np.nonzero(change)[0]]\n    change = change[~np.isnan(change)]\n    change = change[change != -np.inf]\n    change = change[change != np.inf]\n    return np.mean(change)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"percentiles = [1, 5, 10, 20, 25, 30, 40, 50, 60, 70, 75, 80, 90, 95, 99]\nhann_windows = [50, 150, 1500, 15000]\nspans = [300, 3000, 30000, 50000]\nwindows = [10, 50, 100, 500, 1000, 10000]\nborders = list(range(-4000, 4001, 1000))\npeaks = [10, 20, 50, 100]\ncoefs = [1, 5, 10, 50, 100]\nlags = [10, 100, 1000, 10000]\nautocorr_lags = [5, 10, 50, 100, 500, 1000, 5000, 10000]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"W naszym początkowym najprostszym modelu, cechy jak wyznaczyliśmy głównie opierały się na podstawowych miarach tendecji centralnej. W wersji rozszerzonej pozostawiliśmy wcześniejsze cechy(feature) oraz dodaliśmy nowe jak widać poniżej. W skład nowo dodanych cech zaliczają się takie rzeczy jak: dyskretna transformata Fouriera, średnie odchylenie bezwzględne, wartość bezwzględna różnicy średniej, współczynnik zmienności średniej, wartości minimalne, maksymalne, średnie, minimum, maksimum pierwszych, ostatnich pięćdziesięciu tysięcy i dziesięciu tysięcy. Różnice między wartościami minimalnymi i maksymalnymi, współczynnik zmienności średniej dla różnych części zbioru. Dodatkowo użylismy kwantyli, odchylenie standardowe, działania na części rzeczywistej i urojonej szybkiej transformaty Fouriera. Funkcji do przekształcania sygnałów; przekształcenie Hilberta, funkcji Hann, algorytmu STA/LTA. Zastosowaliśmy również wykładniczo ważoną średnią ruchomą. Przeprowadzaliśmy regresje na wybranych podpróbkach. Wykorzystaliśmy funkcję przeprowadzającą tzw. Rolling Window Calculation dla wybranych rozmiarów okna(moving window). \n\nPrzy tworzeniu feature'ów korzystaliśmy z dyskusji na kaggle, oraz z poniższych kerneli oraz własnej wiedzy:\n- https://www.kaggle.com/artgor/even-more-features\n- https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples?scriptVersionId=9803210\n- https://www.kaggle.com/abhishek/quite-a-few-features-1-51\n- https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples\n- https://www.kaggle.com/vettejeep/masters-final-project-model-lb-1-392"},{"metadata":{"trusted":false},"cell_type":"code","source":"def gen_features(x, zero_mean=False):\n    if zero_mean==True:\n        x = x-x.mean()\n    strain = {}\n    strain['mean'] = x.mean()\n    strain['std']=x.std()\n    strain['max']=x.max()\n    strain['kurtosis']=x.kurtosis()\n    strain['skew']=x.skew()\n    zc = np.fft.fft(x)\n    realFFT = np.real(zc)\n    imagFFT = np.imag(zc)\n    strain['min']=x.min()\n    strain['sum']=x.sum()\n    strain['mad']=x.mad()\n    strain['median']=x.median()\n    \n    strain['mean_change_abs'] = np.mean(np.diff(x))\n    strain['mean_change_rate'] = np.mean(np.nonzero((np.diff(x) / x[:-1]))[0])\n    strain['abs_max'] = np.abs(x).max()\n    strain['abs_min'] = np.abs(x).min()\n    \n    strain['avg_first_50000'] = x[:50000].mean()\n    strain['avg_last_50000'] = x[-50000:].mean()\n    strain['avg_first_10000'] = x[:10000].mean()\n    strain['avg_last_10000'] = x[-10000:].mean()\n    \n    strain['min_first_50000'] = x[:50000].min()\n    strain['min_last_50000'] = x[-50000:].min()\n    strain['min_first_10000'] = x[:10000].min()\n    strain['min_last_10000'] = x[-10000:].min()\n    \n    strain['max_first_50000'] = x[:50000].max()\n    strain['max_last_50000'] = x[-50000:].max()\n    strain['max_first_10000'] = x[:10000].max()\n    strain['max_last_10000'] = x[-10000:].max()\n    \n    strain['max_to_min'] = x.max() / np.abs(x.min())\n    strain['max_to_min_diff'] = x.max() - np.abs(x.min())\n    strain['count_big'] = len(x[np.abs(x) > 500])\n           \n    strain['mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n    strain['mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])\n    strain['mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n    strain['mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])\n    \n    strain['q95'] = np.quantile(x, 0.95)\n    strain['q99'] = np.quantile(x, 0.99)\n    strain['q05'] = np.quantile(x, 0.05)\n    strain['q01'] = np.quantile(x, 0.01)\n    \n    strain['abs_q95'] = np.quantile(np.abs(x), 0.95)\n    strain['abs_q99'] = np.quantile(np.abs(x), 0.99)\n    strain['abs_q05'] = np.quantile(np.abs(x), 0.05)\n    strain['abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    for autocorr_lag in autocorr_lags:\n        strain['autocorrelation_' + str(autocorr_lag)] = feature_calculators.autocorrelation(x, autocorr_lag)\n    \n    # percentiles on original and absolute values\n    for p in percentiles:\n        strain['percentile_'+str(p)] = np.percentile(x, p)\n        strain['abs_percentile_'+str(p)] = np.percentile(np.abs(x), p)\n    \n#     strain['trend'] = add_trend_feature(x)\n#     strain['abs_trend'] = add_trend_feature(x, abs_values=True)\n    strain['abs_mean'] = np.abs(x).mean()\n    strain['abs_std'] = np.abs(x).std()\n    \n    strain['quantile_0.95']=np.quantile(x, 0.95)\n    strain['quantile_0.99']=np.quantile(x, 0.99)\n    strain['quantile_0.05']=np.quantile(x, 0.05)\n    strain['realFFT_mean']=realFFT.mean()\n    strain['realFFT_std']=realFFT.std()\n    strain['realFFT_max']=realFFT.max()\n    strain['realFFT_min']=realFFT.min()\n    strain['imagFFT_mean']=imagFFT.mean()\n    strain['imagFFT_std']=realFFT.std()\n    strain['imagFFT_max']=realFFT.max()\n    strain['imaglFFT_min']=realFFT.min()\n    \n    strain['std_first_50000']=x[:50000].std()\n    strain['std_last_50000']=x[-50000:].std()\n    strain['std_first_25000']=x[:25000].std()\n    strain['std_last_25000']=x[-25000:].std()\n    strain['std_first_10000']=x[:10000].std()\n    strain['std_last_10000']=x[-10000:].std()\n    strain['std_first_5000']=x[:5000].std()\n    strain['std_last_5000']=x[-5000:].std()\n        \n    strain['Hilbert_mean'] = np.abs(hilbert(x)).mean()\n    strain['Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n    strain['classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n    strain['classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n    strain['classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n    strain['classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n    #strain['classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean() contains inf and Nan values\n    strain['classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n    #strain['classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean() contains inf and Nan values\n    strain['classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n    strain['Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n    moving_average_700_mean = x.rolling(window=700).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    strain['exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n    strain['exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n    strain['exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n    no_of_std = 3\n    strain['MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n    strain['MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n    \n    strain['iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n    strain['q999'] = np.quantile(x,0.999)\n    strain['q001'] = np.quantile(x,0.001)\n    strain['ave10'] = stats.trim_mean(x, 0.1)\n        \n    for window in windows:\n        x_roll_std = x.rolling(window).std().dropna().values\n        x_roll_mean = x.rolling(window).mean().dropna().values\n        \n        strain['ave_roll_std_' + str(window)] = x_roll_std.mean()\n        strain['std_roll_std_' + str(window)] = x_roll_std.std()\n        strain['max_roll_std_' + str(window)] = x_roll_std.max()\n        strain['min_roll_std_' + str(window)] = x_roll_std.min()\n        strain['q01_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.01)\n        strain['q05_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.05)\n        strain['q95_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.95)\n        strain['q99_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.99)\n        strain['av_change_abs_roll_std_' + str(window)] = np.mean(np.diff(x_roll_std))\n        strain['av_change_rate_roll_std_' + str(window)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n        strain['abs_max_roll_std_' + str(window)] = np.abs(x_roll_std).max()\n        \n        for p in percentiles:\n            strain['percentile_roll_std_' + str(p) + '_window_' + str(window)] = np.percentile(x_roll_std, p)\n            strain['percentile_roll_mean_' + str(p) + '_window_' + str(window)] = np.percentile(x_roll_mean, p)\n        \n        strain['ave_roll_mean_' + str(window)] = x_roll_mean.mean()\n        strain['std_roll_mean_' + str(window)] = x_roll_mean.std()\n        strain['max_roll_mean_' + str(window)] = x_roll_mean.max()\n        strain['min_roll_mean_' + str(window)] = x_roll_mean.min()\n        strain['q01_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.01)\n        strain['q05_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.05)\n        strain['q95_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.95)\n        strain['q99_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.99)\n        strain['av_change_abs_roll_mean_' + str(window)] = np.mean(np.diff(x_roll_mean))\n        strain['av_change_rate_roll_mean_' + str(window)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n        strain['abs_max_roll_mean_' + str(window)] = np.abs(x_roll_mean).max()\n        \n        \n    return pd.Series(strain)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Dalej przeszliśmy do analizy wygenerowanych feature'ów w kontekście tworzenia modelu regresji"},{"metadata":{"trusted":false},"cell_type":"code","source":"del train_df\ngc.collect()\n\ntrain_df = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'), iterator=True, chunksize=150_000, dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\nX_train = pd.DataFrame()\nX_train_zero_mean = pd.DataFrame()\ny_train = pd.Series()\n\nfor df in train_df:\n    features = gen_features(df['acoustic_data'])\n    ch_zero_mean = gen_features(df['acoustic_data'], zero_mean=True)\n    X_train = X_train.append(features, ignore_index=True)\n    X_train_zero_mean = X_train_zero_mean.append(ch_zero_mean, ignore_index=True)\n    y_train = y_train.append(pd.Series(df['time_to_failure'].values[-1]), ignore_index=True)\n\nX_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"W tym miejscu generowaliśmy cechy(features) dla każdego pomiaru danych akustycznych znajdujących się w zbiorze testowym. Dla każdego segmentu danych obliczane były zdefiniowane cechy oraz wyniki przypisywane do struktury X_test."},{"metadata":{"trusted":false},"cell_type":"code","source":"X_test = pd.DataFrame()\nX_test_zero_mean = pd.DataFrame()\n\nfor i, f in enumerate(ld):\n    df = pd.read_csv(os.path.join(TEST_DIR, f))\n    features = gen_features(df['acoustic_data'])\n    ch_zero_mean = gen_features(df['acoustic_data'], zero_mean=True)\n    X_test = X_test.append(features, ignore_index=True)\n    X_test_zero_mean = X_test_zero_mean.append(ch_zero_mean, ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Funkcja odpowiedzialna za rysowanie wykresu przedstawiającego uśrednione dane akustyczne wraz z czasem do trzęsienia dla cech przekazanych w argumentach funkcji. Przy odpowiednim wyborze cech, jesteśmy w stanie po narysowaniu wykresów stwierdzić czy są one przydatne czy nie. W kolejnych krokach będziemy mogli zauważyć co zwraca wywołanie funkcji dla poszczególnych cech, dzięki czemu łatwiej będzie określić ich ważność."},{"metadata":{"trusted":false},"cell_type":"code","source":"def plot_acc_agg_ttf_data(features, title=\"Averaged accoustic data and ttf\"):\n    fig, axes = plt.subplots(3,3, figsize=(30, 18))\n    \n    for i in range(9):\n        plt.title('Averaged accoustic data ({}) and time to failure'.format(features[i]))\n        axes[int(i/3)][i%3].plot(X_train[features[i]], color='r')\n        axes[int(i/3)][i%3].set_xlabel('training samples')\n        axes[int(i/3)][i%3].set_ylabel('acoustic data ({})'.format(features[i]), color='r')\n        plt.legend(['acoustic data ({})'.format(features[i])], loc=(0.01, 0.95))\n        ax2 = axes[int(i/3)][i%3].twinx()\n        ax2.plot(y_train, color='b')\n        ax2.set_ylabel('time to failure', color='b')\n        plt.legend(['time to failure'], loc=(0.01, 0.9))\n        plt.grid(True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Funkcja odpowiedzialna za rysowanie wykresu przedstawiającego histogramy wybranych cech dla zbioru treningowego i testowego, która przydatna jest w stwierdzeniu czy wyniki cech ze zbiorów danych są do siebie zbliżone czy jednak się od siebie różnią. W kolejnych krokach wywołanie funkcji przedstawi nam, że wartości cech zbiorów bardzo od siebie odbiegają, a histogramy prawie nie nachodzą na siebie. Z tego też powodu będziemy starać się zmniejszyć tą rozbieżność i w razie potrzeby przeskalować wartości pomiarów."},{"metadata":{"trusted":false},"cell_type":"code","source":"def plot_distplot_features(features, nlines=4, colors=['green', 'blue'], df1=X_train, df2=X_test):\n    plt.figure()\n    fig, ax = plt.subplots(nlines,2,figsize=(16,4*nlines))\n    for i in range(len(features)):\n        plt.subplot(nlines,2,i+1)\n        plt.hist(df1[features[i]],color=colors[0],bins=50, label='train', alpha=0.5)\n        plt.hist(df2[features[i]],color=colors[1],bins=50, label='test', alpha=0.5)\n        plt.legend()\n        plt.title(features[i])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Jak widać na poniższym wykresie wartości przekazanych cech policzonych dla obu zbiorów odbiegają od siebie. Najbardziej widoczne jest to w przypadku średniej oraz sumy. Na reszczie wykresów rozbieżność jest może troszkę mniejsza lecz wciąż zauważalna."},{"metadata":{"trusted":false},"cell_type":"code","source":"features = ['mean', 'std', 'max', 'min', 'sum', 'mad', 'kurtosis', 'skew']\nplot_distplot_features(features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"features = ['mean', 'std', 'max', 'min', 'sum', 'mad', 'kurtosis', 'skew']\nplot_distplot_features(features, df1=X_train_zero_mean, df2=X_test_zero_mean)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"W przypadku odchylenia standadowego dla różnych części danych sytuacja wygląda tak jak na poniższym wykresie. Wartości dla danych treningowych odbiegają znacznie od danych testowych. Są to różnice wynoszące około 1500-2500."},{"metadata":{"trusted":false},"cell_type":"code","source":"features = ['std_first_50000', 'std_last_50000', 'std_first_25000','std_last_25000', 'std_first_10000','std_last_10000']\nplot_distplot_features(features,3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"features = ['std_first_50000', 'std_last_50000', 'std_first_25000','std_last_25000', 'std_first_10000','std_last_10000']\nplot_distplot_features(features,3, df1=X_train_zero_mean, df2=X_test_zero_mean)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Na poniższym wykresie wzięliśmy losowo 9 featur'ów i następni przedstawiliśmy ich wartości w porównaniu ze zmienną wyjściową 'time_to_failure'\n\nJak widać przy niektórych zmiennych możemy zauważyćpodobieństwo, mianowicie następuje gwałtowny wzrost lub spadek wartości zmiennej razem ze wzrostem wartości zmiennej wyjściowej \"time_to_failure\". "},{"metadata":{"trusted":false},"cell_type":"code","source":"all_features = X_train.columns.values\nnp.random.seed(2019)\nrand_feat_idx = np.random.randint(0, len(all_features), size=9, dtype=np.int32)\nrand_labels = [all_features[x] for x in rand_feat_idx]\n\nplot_acc_agg_ttf_data(rand_labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Postanowiliśmy więc zbadać korelacje między wszystkimi feature'ami a wartością zmiennej wyjściowej\n\nW tym kroku postanowiliśmy wyznaczyć poszczególne korelacje pomiędzy cechami tak, aby zobaczyć, które mają najlepsze wyniki. Pod uwagę braliśmy korelacje o wartości bezwzględnej większej lub równej 0.3. Jak widać w wyniku poniższego kawałku kodu takich korelacji mamy: X. Najlepsze z nich mają wartości powyżej 0.6, to na nich powinniśmy skupić swoją największą uwagę i to one powinny mieć największy wpływ na nasz model."},{"metadata":{"trusted":false},"cell_type":"code","source":"corelations = np.abs(X_train.corrwith(y_train)).sort_values(ascending=False)\ncorelations_df = pd.DataFrame(data=corelations, columns=['corr'])\nprint(\"Number of high corelated values: \",corelations_df[corelations_df['corr']>=0.3]['corr'].count())\n\nhigh_corr = corelations_df[corelations_df['corr']>=0.3]\nprint(high_corr)\nhigh_corr_labels = high_corr.reset_index()['index'].values\nprint(high_corr_labels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"plot_acc_agg_ttf_data(high_corr_labels[:9])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Jak widać na powyższym wykresie nasze przypuszczenia okazały się słuszne i feature'y w których występuję gwałtowny wzrost lub spadek wartości posiadają wysoką korelacje ze zmienną wyjściową"},{"metadata":{},"cell_type":"markdown","source":"Na tym zakończymy proces Analizy i eksploracji danych\n\nW dalszej części projektu przejdziemy do dalszej cześci selekcji oraz przygotowania danych do modelu regresji oraz wyboru samego modelu i jego parametrów"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.7"}},"nbformat":4,"nbformat_minor":1}