目次
今年はどれぐらい降るかねえ
自己紹介
分析の流れ
今後の課題
今年はどれぐらい降るかねえ
私の祖父母がいる地方では、秋が深まるとこんな会話がちらほら。
いわゆる雪深い地域です。
pythonを勉強して約半年ということで
「どれぐらい降るかねえ？」に、いままで勉強してきたことを使って応えてみたいと思います。
自己紹介
2022年春より、Aidemyのデータ分析講座を受講。
勉強してきたことのまとめとして、チューターの先生にサポートしていただきながらこの記事を作りました。
分析の流れ
（１）データの取得、欠損値の処理
（２）データ３種類を1つに結合
（３）それぞれのカラムに前後3日分のデータを追加
（４）線形回帰でモデル実装
（５）決定係数とMSEで評価
（６）LightBGMの実装と評価
はじめに
環境
・Windows11
・google colaboratory
データの引用元
・積雪量データ
この「中越地方>三条」をダウンロードしました。
・気温、降水量、風速データ
2000年~2010年、2010年~2020年の平均気温、降水量、平均風速を選択し、ダウンロード。
（１）データの取得、欠損値の処理
まず必要なモジュールをインポート。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import lightgbm as lgb
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
３つのデータを読み込んで整えます。
まず１つ目。A市の積雪データ。
➀A市の積雪データ
#データを読み込み
df_snow=pd.read_csv("/content/drive/MyDrive/成果物/A市積雪.csv",encoding="shiftjis")
#インデックスを整えます
df_snow = df_snow.set_index('地点13')
#これを新しい変数名へ
snow_new_df = pd.DataFrame()
for i in range(len(df_snow.columns)):
tmp_df = df_snow.iloc[:,i]
tmp_df.index = map( lambda x: str(tmp_df.name)+'年'+str(x),list(tmp_df.index))
tmp_df.columns = '降雪量'
snow_new_df = pd.concat([snow_new_df, tmp_df],axis=0)
snow_new_df.columns = ['降雪量']
出力結果
実は、取り込む前のエクセルデータの時点で一部に「9999」という記述がありました。閏年の関係でしょうか？
これは欠損値として、除きます。
#9999が含まれないようにする
snow_new_df = snow_new_df[snow_new_df['降雪量']<1000]
ですが...閏年以外では2/29が存在しないはずなのに、いくつかの年には2/29に実測値が入っていました。
それも抽出して除きます。
#存在しない日を抽出
except_days = []
for i in snow_new_df.index:
try:
tmp = pd.to_datetime(i, format='%Y年%m月%d日')
except:
except_days.append(i)
これらが、存在しないはずの2/29に実測値が入っていた箇所です。これを除きます。
#除く
snow_new_df = snow_new_df[~snow_new_df.index.isin(except_days)]
#datetime型に整える
snow_new_df.index= pd.to_datetime(snow_new_df.index, format='%Y年%m月%d日')
datetime型に整えたところで、➀積雪データの処理は完了です。
➁A市の降水量・風速・気温データ（2000年～2010年）
読み込んで、欠損値を確認。
#データの読み込み
df_0010=pd.read_csv("/content/drive/MyDrive/成果物/A市00-10 温度・降水量・風速.csv",encoding="shiftjis", skiprows=3)
#欠損値の数をカウント
df_0010.isnull().sum()
#欠損値の確認
df_0010[df_0010['平均気温(℃)'].isnull()]
取り込んだデータ出力
欠損値の数の出力
１はおそらくNaNのこと。
平均気温の４が気になりますね。
#欠損値の確認
df_0010[df_wind_0010['平均気温(℃)'].isnull()]
出力結果
欠損値の箇所を大元のエクセルデータで確認してみると、この11/25-27の気温について欠測（測定がもともと行われていなかった日）があったようでした。
使用するデータ量は足りているので、今回はこの３日間を除きます。
df_0010.dropna(inplace=True)#欠測日をdrop
➂A市の降水量・風速・気温データ（2010年～2020年）
では、同じ要領でもう10年分。
# A市の10-20の風速などのデータを読み込み
df_1020 = pd.read_csv('/content/drive/MyDrive/成果物/v3802/三条市10-20 温度・降水量・風速.csv', encoding='shiftjis', skiprows=3)
出力結果
10年分、とれていますね！
欠損値は...
#欠損値の数をカウント
df_1020.isnull().sum()
１か所。最初のNaNですね。
#NaNをdrop
df_1020.dropna(inplace=True)
dropしておきます。
これでデータが整いました！
（２）データ３種類を1つに結合
これで３つのデータを用意できました！
続いては結合です。
➀A市の降水量・風速・気温データを結合
まずは、ダウンロードの問題で２つに分けて取り込んだデータを１つにします。
#10年分ずつのデータフレームを1つに結合
df_0020=pd.concat([df_0010, df_1020])#結合
df_0020["年月日"]=pd.to_datetime(df_0020["年月日"], format="%Y/%m/%d")#積雪量のデータと形式が合うようにdatetime型へ
➁A市の降水量・風速・気温データと、降雪量データを結合
all_df=pd.merge(snow_new_df, df_0020, left_on=snow_new_df.index, right_on=df_0020["年月日"], how="inner")#結合
all_df=all_df.drop("key_0",axis=1)#インデックス番号は不要なので、消します
all_df= all_df.set_index("年月日")#インデックス番号の代わりに年月日を。
thresh=pd.date_range(start="2000/3/31", end="2020/3/31", freq="12M")
#データを3/31で区切ります。調べるのは毎年12～3月までのデータ。
all_df[(all_df.index>thresh[1]) & (all_df.index<= thresh[2])].sort_index()#2000年から2020年まで、日付を順番に並べ替えます
これで、データを1つにできました！
（３）それぞれのカラムに前後3日分のデータを追加
当日前後の気温や風速なども影響すると思うので、それぞれ3日分のデータを持たせていきます。
勉強したfor文を使って1つずつデータをずらし、前後3日のデータを追加。
all_df=all_df.sort_index()
all_merge_df = pd.DataFrame()
shift_col_name = ['降雪量','平均気温(℃)','平均風速(m/s)']
for i in range(len(thresh)-1):
tmp_df = all_df[(all_df.index > thresh[i]) & (all_df.index <= thresh[i+1])]
for col_name in shift_col_name:#col_name="降雪量"
for num in range(1,4):#num=1
new_col_name = col_name + str(num)#降雪量1
tmp_df[new_col_name] = tmp_df[col_name].shift(num)
all_merge_df = pd.concat([all_merge_df,tmp_df[3:]])
これでデータの前処理が終わりました。
次は、モデル実装です。
（４）線形回帰でモデル実装
今回は、教材でも勉強した線形重回帰モデルを実装します。
lr = LinearRegression()
X = all_merge_df.drop('降雪量',axis=1) # 説明変数（Numpyの配列）
Y = all_merge_df['降雪量'] # 目的変数（Numpyの配列）
train_X, test_X, train_y, test_y = train_test_split(X, Y, random_state=42)# 学習用、テスト用に分割
std_train_X = (train_X - train_X.mean(axis=0)) / train_X.std(axis=0)
std_test_X = (test_X - train_X.mean(axis=0)) / train_X.std(axis=0)#データの標準化
lr.fit(std_train_X, train_y)
# plt.scatter(X, Y, color = 'blue') # 説明変数と目的変数のデータ点の散布図をプロット
pred_y=lr.predict(std_test_X)
mse = mean_squared_error( test_y, pred_y)
print(mse)#MSEの出力
plt.scatter(test_y, pred_y, color = 'red') # 回帰直線をプロット
# plt.title('実測値と予測値のプロット') # 図のタイトル
plt.xlabel('実測値') # x軸のラベル
plt.ylabel('予測値') # y軸のラベル
plt.grid() # グリッド線を表示
plt.show()
グラフの形は、線形回帰のお手本に近めではないかと思います。
雪が降らない日も多いので、0に点が集中していますね。
MSEは31.99...こちらも悪くはない数値です。
（５）決定係数とMSEでモデル評価
モデル評価の方法として、決定係数とMSEを使ってみます。
まずは決定係数。１に近いほど良いとされます。
lr.score(std_test_X, test_y)
だいぶ１に近いですね！
予測と実測値の値が大きく離れていないことを表しています。
MSEの結果は、さきほど出ていたように31.99...でした。
値が小さいほど良いものです。
こちらも、予測と実測値に大きな差がないことを表します。
（６）LightGBMの実装と評価
線形重回帰ともう一つ。
発展として、LightGBMを使って同じ要領で予測してみました。
lgbm = lgb.LGBMRegressor()
X = all_merge_df.drop('降雪量',axis=1) # 説明変数（Numpyの配列）
Y = all_merge_df['降雪量'] # 目的変数（Numpyの配列）
train_X, test_X, train_y, test_y = train_test_split(X, Y, random_state=42)# 学習用、テスト用に分割
lgbm.fit(train_X, train_y)
# plt.scatter(X, Y, color = 'blue') # 説明変数と目的変数のデータ点の散布図をプロット
pred_y=lgbm.predict(test_X)
plt.scatter(test_y, pred_y, color = 'red') # 回帰直線をプロット
mse = mean_squared_error( test_y, pred_y)
print(mse)#MSEの出力
# plt.title('実測値と予測値のプロット') # 図のタイトル
plt.xlabel('実測値') # x軸のラベル
plt.ylabel('予測値') # y軸のラベル
plt.grid() # グリッド線を表示
plt.show()
MSEは41.68...と、線形重回帰よりも10ほど悪くなってしまいました。
lgbm.feature_importances_#特徴重要度を出力
feature_importancesで、どの特徴がどのように予測に影響を与えているのか見てみましょう。
カラムと照らし合わせてみると
1日前の降雪量：503
当日の平均気温：311
この二つの影響力が大きそうですね。
線形重回帰では...？
１日前の降雪量の影響の大きさが目立ちます。
大きな差があって、次点は2日前の平均気温でした。
どちらも、1日前の降雪量にいちばん影響されることには変わりないですね。
どの特徴が予測へ影響するのか、という点では少しモデルの違いが出たということでしょうか。
備考
今回のデータはA市のなかでも雪が比較的少ない地点のものでした。そのため、地元のS地区ではこの予測よりも数値が大きくなることが予想できます。
S地区ピンポイントでの予測も考えたのですが、S地区内での観測データが無く...
A市の気温などを少しづつ下げて、S地区バージョンになおす？とも考えましたが
不確定要素が増えすぎてしまうため断念。
データ分析・予測には何よりデータそのものが大事なのだと、改めて実感しました。
今後の課題
LightGBMについて
こちらは線形回帰よりも形が良くない印象を受けました。
今回は実装にあたりパラメーターの値はデフォルトのまま。ここをうまく調整できたら線形回帰よりも精度が上がるのではないかと思います。
また、雪の降る条件としては上空の温度や、地上の湿度も影響力が大きいです。
今回はデータがなくこの２点は考慮されていないですが、これらが入ればもっと正確な予測ができるのではないかと思います。
データについて
今回使用したデータ自体に0が多かったことについてです。
そもそも雪が降らない日が多ければ実測値、予測値ともに0が多くなります。
そのため決定係数もかなり良い数値が出ていると考えられます。
累積の積雪量予測や、1か月ごとの積雪量予測であれば0がでることはないので、そういった意味では1か月ごと予測をしてみるのも良いのかなと思います。
１日ごとの予測について
今回は１日ごとに予測をしています。ということは、常に前日の実測値から明日の予測をしている。。
その日の実測値から予測できるのは、明日の数値のみ。
今後できるとしたら、今回予測で使った「実測値」を「自分で予測した値」に書き換えてみることでしょうか。
予測した明日のデータを基にあさってを予測➡あさっての予測を基にしあさって...
現時点でのモデル評価指数は良いので、きっとこのような予測でも大きく外れないのではと思います。
さいごに
「どれぐらい降るかねえ？」に明確には答えられないけれど、
python学習のアウトプットという、私にとっては大いに意味のある分析・予測になりました。
今後も、課題点に自分なりに答えをだしていきたいです。
最後に、ここまでサポートして下さったAidemyの先生方に改めて。ありがとうございました！