case-kの備忘録

日々の備忘録です。データ分析とか基盤系に興味あります。

地域の気温をPythonで重回帰分析して予測してみた

今回は地域と時間を説明変数に重回帰分析で地域ごとの気温を重回帰分析で予測します。

本記事の目的

本記事は以下を目的としています。
・重回帰分析の理解
・実際に重回帰分析で問題を解き理解を深める

重回帰分析とは

重回帰分析は1つの目的変数を複数の説明変数で予測する手法です。
例えば、身長から体重を予測するのが単回帰分析で、身長と腹囲と胸囲から体重を予測するのが重回帰分析となります。

実施手順

実際に重回帰分析で問題をといてみたいと思います。
1. サンプルデータの作成
2. 標本分布の確認
3. モデルの学習と評価
4. 交差検証で評価

1. サンプルデータの作成

気象庁の地域・時間ごとの天気をDLし、その後地域には緯度経度を付与します。

# サンプルデータはこちらからDLしてください。
[https://github.com/case-k-git/python/tree/master/%E4%BB%AE%E8%AA%AC%E6%A4%9C%E5%AE%9A/data]

# ライブラリ
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# create sample data
data = pd.DataFrame(columns=columns)
data1 =  pd.read_csv('./data/hokaido_hiro_20180629.csv')
data2 =  pd.read_csv('./data/hokaido_asahikawa_20180629.csv')
data3 =  pd.read_csv('./data/hokaido_kamikawa_20180629.csv')
data4 =  pd.read_csv('./data/tokyo_20180629.csv')
data5 =  pd.read_csv('./data/tokyo_nerima_20180629.csv')
data6 =  pd.read_csv('./data/tokyo_oume_20180629.csv')
data = data.append([data1,data2,data3,data4,data5,data6]) 
data = data[['hour', 'area', 'temperature']]
print(data.head())

# out put
   hour area  temperature
0   1.0   広尾         17.7
1   2.0   広尾         16.2
2   3.0   広尾         16.3
3   4.0   広尾         15.7
4   5.0   広尾         14.5

DLしたデータの地域を調べます。

# unique area data
keys = data['area'][data['area'].duplicated()==False]
print(keys)

# out put 
0    広尾
0    旭川
0    上川
0    東京
0    練馬
0    青梅

緯度経度を付与します。

# area data
columns = ['area','lat', 'long']
values = np.array([['広尾',42.285838,143.311496],
          ['旭川',43.770636,142.364819],
          ['上川',43.84718,142.770431],
          ['東京',35.709026,139.731993],
          ['練馬',35.735623,139.651658],
          ['青梅',35.787996,139.27583],
         ])
area_data = pd.DataFrame(columns=columns,data=values)
print(area_data)

# out put 
  area        lat        long
0   広尾  42.285838  143.311496
1   旭川  43.770636  142.364819
2   上川   43.84718  142.770431
3   東京  35.709026  139.731993
4   練馬  35.735623  139.651658
5   青梅  35.787996   139.27583

areaをキーに緯度経度を気象庁のデータに付与します。

# join area data
area_data = pd.merge(data,area_data, on='area')
print(area_data.head())

# out put 
   hour area  temperature        lat        long
0   1.0   広尾         17.7  42.285838  143.311496
1   2.0   広尾         16.2  42.285838  143.311496
2   3.0   広尾         16.3  42.285838  143.311496
3   4.0   広尾         15.7  42.285838  143.311496
4   5.0   広尾         14.5  42.285838  143.311496
2. 標本分布の確認
# sample distribution
X = area_data['lat']
y = area_data['temperature']
plt.scatter(X, y)
plt.xlabel('lat')
plt.ylabel('temperature')

f:id:casekblog:20180701192006p:plain:w200
エリアによって気温が異なることが確認できます。

3. モデルの学習と評価
columns = ['hour','lat','long']
lr = LinearRegression()
X = area_data[columns]
y = area_data['temperature']
X_train, X_test, y_train, y_test = train_test_split(X,y,random_state=0)
lr.fit(X_train,y_train)
print("test set accuracy:{}".format(lr.score(X_test,y_test)))

# out put
test set accuracy:0.831282607327257

評価精度は83%と非常に高いことが確認できます。
しかし、この数値はサンプルデータを単純に2分割したものです。
交差検証を行いデータの分割を何度も繰り返し行い複数のモデルを訓練する回帰問題のためk分割交差検証を行います。
case-k.hatenablog.com

4. 交差検証で評価
# 交差検証
best_score = 0
for k in range(len(columns)):
    X = area_data[columns].iloc[:,0:k+1]
    mean_score = cross_val_score(lr, X,y,cv=3).mean()
    print("cross validation score:{}, X:{}".format(mean_score,X.keys()))
    if mean_score > best_score:
        best_score = mean_score
print("best score:{}".format(best_score))

# out put 
cross validation score:-3.5037803947710913, X:Index(['hour'], dtype='object')
cross validation score:-0.0724441985185335, X:Index(['hour', 'lat'], dtype='object')
cross validation score:0.04954983428403573, X:Index(['hour', 'lat', 'long'], dtype='object')
best score:0.04954983428403573

交差検証を行った結果、train_test_split関数を活用した時に比べ評価精度が4%と非常に低くなりました。
このモデルは訓練に用いられた特定の分割に強く依存していることを示唆してますが、
今回は単純にデータセットのサイズがちいさすぎるからだと思わます。


以上となります。