Будет ли завтра дождь? Прогнозирование с использованием 5 алгоритмов машинного обучения

В этой статье я расскажу о своем проекте по прогнозированию дождя на следующий день. Я написал 5 алгоритмов (логистическая регрессия, машина опорных векторов, дерево решений, случайный лес, повышение градиента) и сравнил их производительность. Какая из них лучше всего подходит для этой задачи? Продолжайте читать, чтобы узнать.

Источник данных

Первоисточником является Бюро метеорологии правительства Австралии. Файл .csv, использованный в этом анализе, был подготовлен преподавателями курса IBM Машинное обучение с Python. Спасибо Джозефу Сантарканджело, доктору философии. и Светлана Крамар. Файл содержит наблюдения за погодными показателями за каждый день с 2008 по 2017 год в Сиднее, Новый Южный Уэльс, Австралия.

В следующей таблице показаны имена столбцов и их описания. RainTomorrow — наша целевая переменная для классификации.

Примечание. Мой проект написан в Jupyter Notebook.

Импортировать библиотеки

import pandas as pd
import numpy as np

from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV

from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import GradientBoostingClassifier

from sklearn.metrics import classification_report, confusion_matrix, ConfusionMatrixDisplay

import seaborn as sns
import matplotlib.pyplot as plt

Загрузите и поймите набор данных

url = 'https://raw.githubusercontent.com/marvin-rubia/Next-Day-Rain-Prediction/main/Weather_Data.csv'
df = pd.read_csv(url)

df.info()

Результат предыдущего кода:

Поскольку все столбцы имеют одинаковое количество ненулевых значений, это означает, что у нас нет пропущенных значений ни для одного столбца.

Давайте проверим некоторые данные из нашего набора данных, используя df.tail(10):

Предварительная обработка данных

Нам нужно преобразовать строковые значения в числовые фиктивные эквиваленты, потому что нам нужно, чтобы наши функции (или независимые переменные) были нормализованы, прежде чем мы введем их в каждый алгоритм. Чтобы сократить мой блог, я не буду показывать коды предварительной обработки. Но вы можете просмотреть указанные коды в моем блокноте Jupyter, который загружен в мой репозиторий GitHub.

Краткое изложение моих данных таково.

  1. Категориальные «Да» и «Нет» заменены на двоичные 1 и 0 для столбцов RainToday и RainTomorrow.
  2. Столбец Дата преобразован в тип даты.
  3. Добавлен новый столбец под названием Месяц.
    Примечание. Зачем добавлять столбец месяца? Потому что климат зависит от месяцев. Некоторые месяцы более восприимчивы к осадкам, чем другие месяцы.
  4. Столбец Месяц (со значениями от 1 до 12) преобразован в циклические эквиваленты с использованием тригонометрического синуса и косинуса. В результате этого процесса были добавлены два новых столбца: Month_sin и Month_cos.
  5. Преобразовали направления ветра в фиктивные переменные (т.е. горячее кодирование).

После всех этих процессов df.dtypes выдаст следующий результат.

За исключением столбца Дата, все переменные теперь являются числовыми, и именно этот тип наши алгоритмы ожидают в качестве входных данных.

Выбор функций

Осадки вызваны многими факторами окружающей среды. Однако для эффективности, интерпретируемости и, надеюсь, для лучшего прогнозирования я собираюсь использовать свое суждение о науке о дожде.

Чем дальше в будущее мы собираемся предсказывать, тем ниже вероятность того, что мы окажемся точными. Другими словами, если сейчас идет сильный дождь, вполне вероятно, что дождь продолжится еще 10 минут спустя. Но трудно сказать, поможет ли этот пример предсказать, продолжится ли дождь через 7 дней.

Аналогичным образом, даже если дождь еще не, наша уверенность в использовании причинных факторов для прогнозирования когда пойдет дождь уменьшается по мере увеличения времени между моментом наблюдения и моментом прогнозирования. Например, если прямо над моим домом находится нимбовое облако, я могу предсказать, что сегодня пойдет дождь. Но это нимбовое облако не поможет мне предсказать, пойдет ли дождь через 3 дня в моем месте.

Имея это в виду, если мы хотим предсказать, будет ли завтра дождь, наши причинные факторы должны быть ближе к завтрашнему дню. Итак, я собираюсь выбрать наблюдения объектов в 15:00 и исключить их версии в 9 утра.

Кроме того, чтобы повысить эффективность нашего моделирования (избегая ненужных размеров), мы можем исключить WindGustDir и WindGustSpeed, поскольку у нас уже есть WindDir3pm и Скорость ветра15:00. Первая пара описывает направление и скорость самого сильного порыва ветра за день, а вторая пара описывает усредненное направление и скорость ветра до 15:00.

Обратите внимание, что, исходя из этих рассуждений, способность нашей модели предсказывать дождь на следующий день частично зависит от наблюдений в 15:00 предыдущего дня. Наша цель — не построить модель, которая сможет предсказать дальнейшее будущее, а построить модель, которая сможет предсказать возможность дождя завтра, используя любые измерения, которые можно провести во время время работы (в метеообсерватории) с предыдущего дня.

Итак, теперь у нас есть уменьшенный фрейм данных:

# Create the reduced dataframe (selecting only the columns we need for analysis)

# NOTE: We remove Date and Month columns since we already have the cyclical encoding via Month_sin and Month_cos
# NOTE: We rearrange the order of some columns and move the target column (RainTomorrow) to the last column

df_reduced = df[['Month_sin', 'Month_cos', 'MinTemp', 'MaxTemp', 
                 'Rainfall', 'Evaporation', 'Sunshine', 'Cloud3pm', 
                 'Humidity3pm', 'Pressure3pm', 'Temp3pm', 'WindSpeed3pm', 
                 'WindDir3pm_E', 'WindDir3pm_ENE', 'WindDir3pm_ESE', 'WindDir3pm_N', 
                 'WindDir3pm_NE', 'WindDir3pm_NNE', 'WindDir3pm_NNW', 'WindDir3pm_NW', 
                 'WindDir3pm_S', 'WindDir3pm_SE', 'WindDir3pm_SSE', 'WindDir3pm_SSW', 
                 'WindDir3pm_SW', 'WindDir3pm_W', 'WindDir3pm_WNW', 'WindDir3pm_WSW', 
                 'RainToday', 'RainTomorrow']] 

df_reduced.tail(10)

Вывод предыдущих кодов:

Что еще удалить?
Мы можем удалить переменные, которые сильно коррелируют с другими переменными (чтобы избежать мультиколлинеарности). Для наших алгоритмов важно избежать переобучения из-за избыточной информации, и это верно не только для регрессионного моделирования, но и для классификации.

Для этого я собираюсь закодировать корреляционные матрицы для парных переменных с коэффициентами корреляции r ≥ +0,60 и r ≤ -0,60, а затем удалить по одному компоненту из каждой пары.

# Show variables that are highly correlated with r >= +0.6

corr_matrix = df_reduced.corr()

condition =  corr_matrix >= 0.6

plt.figure(figsize=(16, 14))
sns.heatmap(corr_matrix[condition], annot=True, cmap='Blues')
plt.show()
# Show variables that are highly correlated with r <= -0.6

corr_matrix = df_reduced.corr()

condition =  corr_matrix <= -0.6

plt.figure(figsize=(16, 14))
sns.heatmap(corr_matrix[condition], annot=True, cmap='Blues')
plt.show()

Если вы запустите предыдущие коды, вы легко увидите, какие столбцы сильно коррелируют. В результате этого процесса мы пришли к выводу, что можем удалить MinTemp, MaxTemp, Evaporation и Sunshine.

Теперь у нас есть окончательный фрейм данных для анализа:

df_final = df_reduced.drop(['MinTemp', 'MaxTemp', 'Evaporation', 'Sunshine'], axis=1)

Со следующими характеристиками:

Последняя проверка: есть ли у нас функция с равномерным распространением?

# Show histogram for each feature

df_final.drop('RainTomorrow', axis=1).hist(figsize=(24,24))
plt.show()

Не существует функций с равномерным распределением, поэтому удалять дальше нечего. (Равномерное распределение указывает на то, что независимая переменная имеет постоянную вероятность для любого своего значения, что не поможет нашей модели определить, является ли переменная ответа это 1 или 0.)

Зачем нам вообще нужны направления в наших функциях?
Синди — столица Нового Южного Уэльса и самый густонаселенный город в Австралии. С восточной стороны он окружен южной частью Тихого океана.

В глобальном масштабе пассаты влияют на времена года на Земле. На местном уровне на направление ветра влияют муссоны, которые влияют на количество осадков в определенные месяцы. Итак, да, нам нужны характеристики по направлениям. Особая закономерность в направлении ветра в Сиднее может помочь объяснить выпадение осадков.

Сбалансированный или несбалансированный набор данных?
Давайте проверим записи положительных и отрицательных классов для нашей целевой переменной (RainTomorrow).

# Count the target's distinct values

df['RainTomorrow'].value_counts()

Вывод предыдущего кода показывает, что имеется 2422 записи для отрицательного класса (0 = завтра дождя не будет) и 849 записей для положительного класса (1 = дождь завтра). Только 35 % записей нашей цели относятся к положительному классу, а это означает, что наш набор данных несбалансирован.

Отдельные функции и цель

features = df_reduced.drop('RainTomorrow', axis=1)
y = df['RainTomorrow']

Нормализация функций

Прежде чем приступить к моделированию, нам необходимо нормализовать масштабы выбранных функций, чтобы наши алгоритмы обрабатывали каждый из них одинаково.

scaler = StandardScaler()

normalized_array = scaler.fit_transform(features)
features_scaled = pd.DataFrame(normalized_array, columns=features.columns)

features_scaled.head()

Результат предыдущих кодов:

Создавайте обучающие и тестовые наборы

X_train, X_test, y_train, y_test = train_test_split(features_scaled, y, test_size=0.20, random_state=8)

print(len(X_train), ':', len(y_train))
print(len(X_test), ':', len(y_test))

Если вы запустите предыдущие коды, вы увидите, что X_train и y_train имеют по 2616 записей (как и должно быть), а также X_test и . y_test оба имеют 655 записей (как и должно быть).

Моделирование машинного обучения

Завтра мы будем использовать 5 моделей для классификации дождя (или отсутствия дождя) с использованием логистической регрессии, машины опорных векторов, дерева решений, случайного леса и повышения градиента.

1. Логистическая регрессия

Обучение модели с помощью обучающего набора:

# Create our Logistic Regression model

log_reg = LogisticRegression(solver='liblinear')

# Find the best C parameter using GridSearchCV

param_grid_LR = {'C': [0.001, 0.01, 0.1, 1, 10, 100]}

grid_search_LR = GridSearchCV(log_reg, param_grid_LR, cv=5)

grid_search_LR.fit(X_train, y_train)

# Print the best parameter and accuracy score

print("Best C:", grid_search_LR.best_params_['C'])
print("Best score:", grid_search_LR.best_score_)

# Assign the best Logistic Regression model

log_reg_final = grid_search_LR.best_estimator_
log_reg_final

Оценка нашей модели с помощью тестового набора:

# Evaluate our model

y_pred_LR = log_reg_final.predict(X_test)

report = classification_report(y_test, y_pred_LR)

print(report)

Предыдущие процессы моделирования (т. е. обучение наилучшему оценщику и последующая оценка с помощью тестового набора) просто повторятся для остальных моделей.

2. Машина опорных векторов

# Create our SVM model

svm = SVC()

# Find the best C parameter using GridSearchCV

param_grid_svm = {'C': [0.001, 0.01, 0.1, 1, 10, 100]}

grid_search_svm = GridSearchCV(svm, param_grid_svm, cv=5)

grid_search_svm.fit(X_train, y_train)

# Print the best parameter and accuracy score

print('Best C:', grid_search_svm.best_params_['C'])
print('Best score:', grid_search_svm.best_score_)

# Assign the best SVM model

svm_final = grid_search_svm.best_estimator_
svm_final

# Evaluate our model

y_pred_svm = svm_final.predict(X_test)

report = classification_report(y_test, y_pred_svm)

print(report)

3. Классификатор дерева решений

# Create our DTC model

dtc = DecisionTreeClassifier()

# Find the best max_depth parameter using GridSearchCV

param_grid_dtc = {'max_depth': list(range(1, 16))}

grid_search_dtc = GridSearchCV(dtc, param_grid_dtc, cv=5)

grid_search_dtc.fit(X_train, y_train)

# Print the best parameter and accuracy score

print('Best max_depth:', grid_search_dtc.best_params_['max_depth'])
print('Best score:', grid_search_dtc.best_score_)

# Assign our best DTC model

dtc_final = grid_search_dtc.best_estimator_
dtc_final

# Evaluate our model

y_pred_dtc = dtc_final.predict(X_test)

report = classification_report(y_test, y_pred_dtc)

print(report)

4. Классификатор случайного леса

# Create our RFC model

rfc = RandomForestClassifier(random_state=8)

# Find the best max_features parameter using GridSearchCV

param_grid_rfc = {'max_features': list(range(1,26))}

grid_search_rfc = GridSearchCV(rfc, param_grid_rfc, cv=5)

grid_search_rfc.fit(X_train, y_train)

# Print the best parameter and accuracy score

print('Best max_features:', grid_search_rfc.best_params_['max_features'])
print('Best score:', grid_search_rfc.best_score_)

# Assign our best RFC model

rfc_final = grid_search_rfc.best_estimator_
rfc_final

Из-за загрузочной природы алгоритма случайного леса вы можете получить немного значение лучшего результата (и, следовательно, значение лучшего max_featuresпараметр). Случайный лес использует случайные загрузочные выборки (с заменой) из исходного набора данных для построения отдельных деревьев решений.

Для каждого дерева при каждом запуске используется немного разное подмножество, поэтому производительность немного варьируется. Самый простой способ противостоять этому и обеспечить согласованность — установить для параметра random_state определенное значение. Точная стоимость не имеет значения. (В предыдущих кодах я выбрал 8, поскольку родился 8 сентября на Филиппинах.)

# Evaluate our model

y_pred_rfc = rfc_final.predict(X_test)

report = classification_report(y_test, y_pred_rfc)

print(report)

5. Классификатор повышения градиента

# Create our GBC model

gbc = GradientBoostingClassifier()

# Find the best learning_rate parameter using GridSearchCV

param_grid_gbc = {'learning_rate': [0.001, 0.01, 0.1, 0.2]}

grid_search_gbc = GridSearchCV(gbc, param_grid_gbc, cv=5)

grid_search_gbc.fit(X_train, y_train)

# Print the best parameter and accuracy score

print('Best learning_rate:', grid_search_gbc.best_params_['learning_rate'])
print('Best score:', grid_search_gbc.best_score_)

# Assign our best GBC model

gbc_final = grid_search_gbc.best_estimator_
gbc_final

# Evaluate our model

y_pred_gbc = gbc_final.predict(X_test)

report = classification_report(y_test, y_pred_gbc)

print(report)

Сравнение моделей

Поскольку наш набор данных несбалансирован (только 35 % записей указывают на дождь завтра), давайте оценивать наши модели, используя взвешенный показатель f1.

В случае равенства мы сравниваем их взвешенные оценки. Зачем вспоминать (а не точность)? Потому что мы хотим минимизировать ложноотрицательные результаты. Другими словами, если мы прогнозируем, что завтра не будет дождя, люди не возьмут с собой зонтики, поэтому нам лучше быть более уверенными в этом прогнозе.

# Summarize scores for our 5 models

data = {'LogReg': [0.83, 0.84],
        'RF': [0.81, 0.82],
        'GBC': [0.80, 0.81],
        'SVM': [0.79, 0.81],
        'DTC': [0.77, 0.79]
        }

clf_scores = pd.DataFrame(data)
clf_scores.index = ['weighted f1', 'weighted recall']
clf_scores

Нашей лучшей моделью является Логистическая регрессия, за которой следует Случайный лес. Между тем, модель Дерева решений оказалась самой слабой для нашего набора данных.

# Create the confusion matrix for our Logistic Regression model

cm = confusion_matrix(y_test, y_pred_LR)

cm_display = ConfusionMatrixDisplay(cm)
cm_display.plot()
plt.title('Confusion Matrix for Logistic Regression')
plt.show()

Заключение

Из пяти моделей, проанализированных для исторических данных о погоде в Сиднее с 2008 по 2017 год, нашей лучшей моделью для прогнозирования дождя на следующий день является Логистическая регрессия (C=0,1, Solver='liblinear', остальная часть параметры по умолчанию). Его взвешенный показатель f1 составляет 0,83. Дерево решений (max_length=4, остальные параметры по умолчанию) выполнило самый слабый прогноз с оценкой 0,77.

Последний комментарий: Как это соотносится с прогнозами погоды по телевидению?

Для справки: нынешние обсерватории имеют точность прогнозов 90 % не только на завтра, но и на 5-дневный прогноз. Это означает, что их прогноз на следующий день может иметь точность более 90%.

Честно говоря, эти прогнозы основаны на спутниковых изображениях, отслеживаемых в реальном времени. Напротив, в нашем наборе данных погодные показатели записаны только до 15:00 предыдущего дня.

Тем не менее, взвешенная оценка f1 0,83, взвешенная полнота 0,84 и точность 0,84 для нашей модели Логистическая регрессия являются хорошими показателями для прогнозирования завтрашнего дождя. Житель Сиднея, который регулярно доверяет нашим прогнозам о том, что на следующий день не будет дождя (и не возьмет с собой зонтик), не промокнет под дождем и будет проклинать нас примерно в 80% случаев.

На простом английском языке

Спасибо, что вы являетесь частью нашего сообщества! Прежде чем уйти: