目次

はじめに
 目的
 開発環境
 使用データ
 予測モデル作成
 予測結果
 終わりに
 参考文献

はじめに

現在エンジニア転職に向け日夜プログラミングの勉強に励んでいるmayamiyaと申します。
今回は私の通っているプログラミングスクール''Aidemy''さんの最終課題として、本記事を投稿します。
至らぬ点が多々あると思いますが、ご指摘や感想等あればコメントしていただければ嬉しいです。

目的

今回は毎年私が苦しんでいる「花粉」について予測できるモデルを作り生活に役立てようと考えます。

開発環境

Mac_Apple M1
Google Colaboratory
Python3.9.16

使用データ

環境省花粉観測システム
・花粉飛散量（さいたま）※2012年から2020年のデータを使用
国土交通省　気象庁
・気温、風速、降水量 (さいたま)※2012年から2020年のデータを使用

予測モデル作成

ライブラリのimport

今回利用するライブラリをimportします。

# import
import glob
import lightgbm as lgb
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import random
import seaborn as sns
import tensorflow as tf

from sklearn.impute import KNNImputer
from sklearn.metrics import mean_absolute_error
from sklearn.model_selection import KFold
from tensorflow.keras.layers import Dense, Activation
from tensorflow.keras.models import Sequential
from tensorflow.keras.utils import to_categorical
from tensorflow.keras.utils import plot_model

データの取得

あらかじめgoogleドライブに保存しておいたデータを取得します。今回はデータ保存の段階で訓練データ(2012~2013)とテストデータ(2020)を分けています。

# データの取得

# 花粉飛散量訓練データ(さいたまのデータだけにする)
df_train = pd.read_csv('/content/drive/MyDrive/pollen_csv/2012_2019_pollen.csv')
df_train = df_train.loc[:, ['年', '月', '日', '時', 'さいたま']]

# 花粉飛散量テストデータ(さいたまのデータだけにする)
df_test = pd.read_csv('/content/drive/MyDrive/pollen_csv/2020_pollen.csv')
df_test = df_test.loc[:, ['年', '月', '日', '時', 'さいたま']]

# 特徴量(気象データ)の訓練データ
# 2012~2019年のweather_dataを合わせる
weather_data = os.path.join('/content/drive/MyDrive/pollen_csv', '*_weather.csv')

weather_datas = sorted(glob.glob(weather_data))
weather_datas = weather_datas[0:-1]

df_weather_train = pd.DataFrame()
for data in weather_datas:
  csv = pd.read_csv(data, encoding='shift_jis',
                         header=2)
  #必要なデータのみを取り出しておく
  csv = csv.drop(csv.index[[0, 1]])
  csv = csv.drop(csv.columns[[5, 6, 8, 9, 10, 11, 13, 14]], axis=1)
  df_weather_train = pd.concat([df_weather_train, csv], ignore_index=True)

# 特徴量(気象データ)のテストデータ
df_weather_test = pd.read_csv('/content/drive/MyDrive/pollen_csv/2020_weather.csv', encoding='shift_jis',
                         header=2)

df_weather_test = df_weather_test.drop(df_weather_test.index[[0, 1]])

df_weather_test = df_weather_test.drop(df_weather_test.columns[[5, 6, 8, 9, 10, 11, 13, 14]], axis=1)

データの前処理

データを確認・利用しやすいように一度処理をします。まずは取得したそれぞれのデータを処理しやすいように結合しカラム'Test_Flag'に0(訓練データ)、1(テストデータ)を挿入する。

# 花粉飛散量の訓練データとテストデータを結合しall_dfを作成
all_df = pd.concat([df_train, df_test],axis=0).reset_index(drop=True)

# 訓練データとテストデータを判定するためのカラムを挿入
all_df['Test_Flag'] = 0
all_df.loc[df_train.shape[0]: , 'Test_Flag'] = 1

all_df

# 特徴量の訓練データとテストデータを結合しall_weather_dfを作成
all_weather_df = pd.concat([df_weather_train, df_weather_test],axis=0).reset_index(drop=True)

# 訓練データとテストデータを判定するためのカラムを挿入
all_weather_df['Test_Flag'] = 0
all_weather_df.loc[df_weather_train.shape[0]: , 'Test_Flag'] = 1

all_weather_df

　all_dfの表を見るとカラム名が「さいたま」になっているので「花粉飛散量」に変更する。

# カラム名「さいたま」を「花粉飛散量」に変更する
all_df = all_df.rename(columns={'さいたま': '花粉飛散量'})

all_df.columns

データの確認

データの中身を統計量から確認する。

# 統計量
print('---統計量---')
print(all_df.describe())

print(all_weather_df.describe())

統計量からわかったことは

花粉飛散量、気温、風速、降水量に欠損値がある
花粉飛散量の最小値・最大値ともに外れ値がある
である。

外れ値と欠損値の処理

確認したデータにあった外れ値と欠損値を処理する。今回は花粉飛散量：-9998.000000, -9997.000000, -9996.000000, 5694.000000を外れ値として処理する。

# 花粉飛散量の外れ値を欠損値に変換
all_df = all_df.replace([-9998.000000, -9997.000000, 
                         -9996.000000, 5694.000000], np.nan)

続いて欠損値の処理を行う。今回は.interpolate()を使って欠損値を補完する。

# 欠損値を補完
all_df = all_df.interpolate()
all_weather_df = all_weather_df.interpolate()
print(all_df.describe())
print(all_weather_df.describe())

データを確認すると欠損値がなくなっていることがわかる。花粉飛散量の最大値が3636と高い数値であるが、3月の花粉飛散量は総じて高く3000~の数値もいくつか出ていたため今回はこのまま予測を立ててみることにする。

特徴量の追加

今の特徴量は種類が少ないためいくつか特徴量を追加する。今回追加する特徴量は前日花粉飛散量、直近24時間および1週間の平均気温・平均降水量・平均風速、1週間の平均花粉飛散量を追加する。
はじめに前日花粉飛散量を追加する。

# 前日花粉飛散量関数
def pre_data(y):
  pre_df = all_df[all_df['年'] == y]['花粉飛散量']
  pre_df = pre_df.drop(pre_df.index[-24:])
  
  # 1月31日データはないがほぼ花粉はないと考え毎年のデータの始まりである2月1日の前日量は０にする
  list = [0 for i in range(24)] 
  
  pre_0 = pd.Series(data=list)
  pre_d = pd.concat([pre_0, pre_df], ignore_index=True)
  return pre_d

# pre_dataを使って前日花粉飛散量を追加
pre_pollen = pd.Series()
year_list = [2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020]
for i in year_list:
  pre_pollen = pd.concat([pre_pollen, pre_data(i)], ignore_index=True)

all_weather_df['前日花粉飛散量'] = pre_pollen

続いて直近24時間および1週間の平均気温・平均降水量・平均風速を追加する。また、モデル(LightGBM)実施の際にエラーを起こさないようカラム名から「()」を削除しておく。

# 「()」の削除
all_weather_df = all_weather_df.rename(columns={'気温(℃)': '気温', '降水量(mm)':'降水量', '風速(m/s)':'風速'})

# 直近24時間の平均気温
all_weather_df["24h平均気温"]=all_weather_df["気温"].rolling(24).mean()

# 直近24時間の平均降水量
all_weather_df["24h平均降水量"] =all_weather_df["降水量"].rolling(24).mean()

# 直近24時間の平均風速
all_weather_df['24h平均風速'] = all_weather_df['風速'].rolling(24).mean()

# 直近1週間の平均気温
all_weather_df["1w_平均気温"]=all_weather_df["気温"].rolling(168).mean()

# 直近1週間の平均降水量
all_weather_df["1w_平均降水量"] =all_weather_df["降水量"].rolling(168).mean()

# 直近1週間の平均風速
all_weather_df["1w_平均風速"] =all_weather_df["風速"].rolling(168).mean()

# 直近1週間の平均花粉飛散量
all_weather_df["1w_平均花粉飛散量"] = all_df["花粉飛散量"].rolling(168).mean()

.rollingを利用すると頭の数値が欠損値になるため、欠損値の処理を行う。前述した.interpolate()は前後の値を参照して欠損値を補完するため今回は使えない。そのためKNNImputer()を用いて欠損値を補完する。

from sklearn.impute import KNNImputer
imputer = KNNImputer(n_neighbors=2)
all_weather_df = pd.DataFrame(imputer.fit_transform(all_weather_df), columns = all_weather_df.columns)

all_weather_df

さらに、月・日・時も特長量として使用したいため、@shimopino(日鉄ソリューションズ株式会社)の記事を利用させていただき、月・日。時を循環性を持つ特長量に変換する。

# 月、日、時を周期性をsin,cosで表す
def encode(df, col):
    df[col + '_cos'] = np.cos(2 * np.pi * df[col] / df[col].max())
    df[col + '_sin'] = np.sin(2 * np.pi * df[col] / df[col].max())
    return df

all_weather_df = encode(all_weather_df, '月')
all_weather_df = encode(all_weather_df, '日')
all_weather_df = encode(all_weather_df, '時')

all_weather_df

LightGBMモデル作成

ようやく待ちに待った予測モデルの作成に移ります。本記事ではモデリングをLightGBMとニューラルネットワーク(MLP)を用いて行います。また、今回はK分割交差検証(KFold)を用いてデータを分割する。
始めにLightGBM(regression_L1)を用いてモデリングをする。

# データを訓練用データと検証用データに分割
train = all_weather_df[all_weather_df['Test_Flag']==0]
test = all_weather_df[all_weather_df['Test_Flag']==1]

train_target = all_df[all_df['Test_Flag']==0]['花粉飛散量']
test_target = all_df[all_df['Test_Flag']==1]['花粉飛散量']

col = ['年', '月', '日', '時', 'Test_Flag']
train = train.drop(col, axis=1)
test = test.drop(col, axis=1)


# K-Foldを定義
cv = KFold(n_splits=3, random_state=0, shuffle=True)

train_acc_list = []
val_acc_list = []
models = []

# LightGBMのパラメータを定義
lgb_params = {
    'objective': 'regression_l1',
    'metric': 'mean_absolute_error',
    'force_row_wise': True,
    'seed': 0,
    'learning_rate': 0.1,
    'min_data_in_leaf': 5
}

# モデルを学習、交差検証法で評価
for i ,(trn_index, val_index) in enumerate(cv.split(train, train_target)):

  print(f'Fold: {i}')
  X_train, X_val = train.iloc[trn_index], train.iloc[val_index]
  y_train, y_val = train_target[trn_index], train_target[val_index]

  # LightGBM用のデータセットを定義
  lgb_train = lgb.Dataset(X_train, y_train)
  lgb_valid = lgb.Dataset(X_val, y_val)

  model = lgb.train(
      params = lgb_params,
      train_set = lgb_train,
      valid_sets = [lgb_train, lgb_valid],
      verbose_eval = 100,
      early_stopping_rounds=10
  )

  y_pred = model.predict(X_train)
  train_acc = mean_absolute_error(y_train, y_pred)
  train_acc_list.append(train_acc)

  y_pred_val = model.predict(X_val)
  val_acc = mean_absolute_error(y_val, y_pred_val)
  val_acc_list.append(val_acc)

  models.append(model)

print('-'*10 + 'Result' + '-'*10)
print(f'Train_acc : {train_acc_list}, Ave : {np.mean(train_acc_list)}')
print(f'Valid_acc : {val_acc_list}, Ave : {np.mean(val_acc_list)}')
print('-'*10 + 'Result' + '-'*10)

# 2020年の花粉飛散量
test_pred = model.predict(test)
test_acc = mean_absolute_error(test_target, test_pred)
print(f'Test_acc : {test_acc}')

# グラフ作成準備
pred = pd.DataFrame(test_pred, columns=['pred'])

test['test_target'] = test_target
test = test.reset_index()
test_result = pd.concat([test, pred], axis=1)

# 予測値及び正解の値（実測値）をグラフに出力
test_result['test_target'].plot(grid=True, color='b')
test_result['pred'].plot(grid=True, color='r')
plt.show()

ニューラルネットワーク(MLP)モデル作成

続いてニューラルネットワーク(MLP)によるモデリングを行います。データの分割は同じくK分割交差検証(KFold)を用います。

# データを訓練用データと検証用データに分割
train = all_weather_df[all_weather_df['Test_Flag']==0]
test = all_weather_df[all_weather_df['Test_Flag']==1]


train_target = all_df[all_df['Test_Flag']==0]['花粉飛散量']
test_target = all_df[all_df['Test_Flag']==1]['花粉飛散量']

col = ['年', '月', '日', '時', 'Test_Flag']
train = train.drop(col, axis=1)
test = test.drop(col, axis=1)

# 乱数を固定
tf.random.set_seed(0)
np.random.seed(0)
random.seed(0)
os.environ["PYTHONHASHSEED"] = "0"

cv = KFold(n_splits=3, random_state=0, shuffle=True)

train_acc_list = []
val_acc_list = []

# fold毎に学習データのインデックスと評価データのインデックスが得られます
for i ,(trn_index, val_index) in enumerate(cv.split(train, train_target)):
    print(f'Fold : {i}')
    # データ全体(Xとy)を学習データと評価データに分割
    X_train ,X_val = train.loc[trn_index], train.loc[val_index]
    y_train ,y_val = train_target[trn_index], train_target[val_index]

    # モデルを定義
    model = Sequential()
    
    model.add(Dense(32, input_dim=17))
    model.add(Activation("relu"))
    model.add(Dense(16))
    model.add(Activation("relu"))
    model.add(Dense(1))

    model.compile(loss='mean_absolute_error', optimizer='adam')
    
    model.fit(X_train, y_train, epochs=20, batch_size=16, verbose=None)

    y_pred = model.predict(X_train)
    train_acc = mean_absolute_error(y_train, y_pred)
    print(train_acc)
    train_acc_list.append(train_acc)
    
    y_pred_val = model.predict(X_val)
    val_acc = mean_absolute_error(
        y_val, y_pred_val
        )
    print(val_acc)
    val_acc_list.append(val_acc)

print('-'*10 + 'Result' +'-'*10)
print(f'Train_acc : {train_acc_list} , Ave : {np.mean(train_acc_list)}')
print(f'Valid_acc : {val_acc_list} , Ave : {np.mean(val_acc_list)}')

# 2020年の花粉飛散量
test_pred = model.predict(test)
test_acc = mean_absolute_error(test_target, test_pred)
print(f'Test_acc : {test_acc}')

# グラフ作成準備
pred = pd.DataFrame(test_pred, columns=['pred'])

test['test_target'] = test_target
test = test.reset_index()
test_result = pd.concat([test, pred], axis=1)

# 予測値及び正解の値（実測値）をグラフに出力
test_result['test_target'].plot(grid=True, color='b')
test_result['pred'].plot(grid=True, color='r')
plt.show()

予測結果

予測結果をまとめると、

＝LightGBMモデル＝
訓練用データの絶対平均誤差（の平均）：19.95
バリデーションデータの絶対平均誤差（の平均）：21.90
テストデータの絶対平均誤差：8.66

＝ニューラルネットワーク(MLP)モデル＝
訓練用データの絶対平均誤差（の平均）：24.09
バリデーションデータの絶対平均誤差（の平均）：24.31
テストデータの絶対平均誤差：7.75

でした。
ニューラルネットワークモデルの方がテストデータを少ない誤差で予測することができました。しかしどちらのモデルも花粉飛散量の大まかな推移や上昇タイミングを読み解くことはできていますが、大きな上昇値を予測することができませんでした。

終わりに

今回のデータ分析では、データ分析の作業はデータの取得・処理がものすごく大変であるということを学ぶことができました。スクーリングの授業ではあらかじめデータが用意されていたので実際に自分で収集・修正・追加していくことに多くの時間を使いました。
ここからのより誤差を少なくしていくための作業は後日行い、少しずつ更新いていこうと思います。今回のデータ収集では見つからなかった気圧や湿度を入れることでより精度の高いモデルを作成できるのではないかと思います。
作業中に起きたエラーや精度を上げる作業ではAidemyの先生方から教わることができ、データ分析について良い復習・アウトプットができました。今後も他のデータ分析を行いながら新たな知識を学んでいきたいと思います。

参考文献

環境省花粉観測システム
 国土交通省　気象庁
 時系列コンペでよく使用されるsin, cosへの変換とは