(Python)初心者のメモ帳
今年の４月からPythonの学習を始めた初心者です。
コロナについて一般人が得られる情報をもとに分析を行いました。
今回は下記のHPから感染者数と関係しそうなデータを下記HPからデータを取得しました。
１．データ取得
・東京都のコロナ感染者数：札幌医科大学HP
(※)今回は「７日後増加人数能力100万人当たりの規格化値」を使いました。
https://web.sapmed.ac.jp/canmol/coronavirus/index.html
・東京における移動傾向：Apple HP
（※）今回は東京における人々の移動量(車なども含む)として使用
https://covid19.apple.com/mobility
・日経平均株価(景況感の情報として)
https://indexes.nikkei.co.jp/nkave/index?type=download
・日経電子版_マーケットのtwitterのネガポジ分析(景況感の情報として)
・気象：気象庁HP
https://www.data.jma.go.jp/gmd/risk/obsdl/index.php
1.1. 感染者数と移動量の相関
メディアの報道によると感染から発症まで２週間程度かかるらしい。
感染者数に対して１週間づつ移動量をずらして決定係数を比較しました(図左)。感染者数の相関係数を調べると２週よりも３,4週ずらした方が相関が高かった。今回はずらしすぎを懸念して3週を採用した。
（補記：移動量はcity,sub_region(23区内外)とdrive,walking,transit(交通手段)の組み合わせで6種類あったが、それたの平均値を用いた）
1.2. 株価
ダウンロードしたcsvは日付に対して番号を振ったダミーファイルとmergeさせた。取引のない日(＝欠損値)は直前の値を参照するffillで埋めました。
(参考)https://qiita.com/0NE_shoT_/items/8db6d909e8b48adcb203
import pandas as pd
#dataset
df_1 = pd.read_csv('dummy.csv',index_col=0)#0列:日付、1列:1,2,3,,,
df_2 = pd.read_csv('nikkei_stock.csv',index_col=0)
#result
result = pd.merge(df_1, df_2, on ='date',how='outer')
result = result.iloc[0:247,1].fillna(method = 'ffill')
result.to_csv('nikkei_stock_ffill.csv')1.3. 日経マーケットのツイート取得とネガポジ分析
景況感が行動及び感染者に影響を与えると考え、日経マーケットのツイートのネガポジ分析を行いました。（※）コードはAidemyの教材の引用です。
・ツイートの取得 (twitter APIが利用申請が必要です。)
import tweepy
import csv
consumer_key = "AAAA" #Consumer Key
consumer_secret = "BBBB" # Consumer Secret
access_key = "CCCC" # Access Token
access_secret = "DDDD" # Accesss Token Secert
auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_key, access_secret)
api = tweepy.API(auth)
# ツイート取得
tweet_data = []
for tweet in tweepy.Cursor(api.user_timeline,screen_name = "@nikkei_market",exclude_replies = True).items():
tweet_data.append([tweet.id,tweet.created_at,tweet.text.replace('\n','')])
# csv出力
with open('tweets_nikkei.csv', 'w',newline='',encoding='utf-8') as f:
writer = csv.writer(f, lineterminator='\n')
writer.writerow(["id","text","created_at"])
writer.writerows(tweet_data)取得したツイートをネガポジ分析
import matplotlib.pyplot as plt
import pandas as pd
import MeCab
import re
m = MeCab.Tagger('')
def get_diclist(text):
parsed = m.parse(text)
lines = parsed.split('\n')
lines = lines[0:-2]
diclist = []
for word in lines:
l = re.split('\t|,',word)
d = {'Surface':l[0], 'POS1':l[1], 'POS2':l[2], 'BaseForm':l[7]}
diclist.append(d)
return(diclist)
import pandas as pd
pn_df = pd.read_csv('pn_ja.csv', encoding='utf-8', names=('Word','Reading','POS', 'PN'))
word_list = list(pn_df['Word'])
pn_list = list(pn_df['PN'])
pn_dict = dict(zip(word_list, pn_list))
import numpy as np
def add_pnvalue(diclist_old, pn_dict):
diclist_new = []
for word in diclist_old:
base = word['BaseForm'] # 個々の辞書から基本形を取得
if base in pn_dict:
pn = float(pn_dict[base])
else:
pn = 'notfound' # その語がPN Tableになかった場合
word['PN'] = pn
diclist_new.append(word)
return(diclist_new)
def get_mean(dictlist):
pn_list = []
for word in dictlist:
pn = word['PN']
if pn!='notfound':
pn_list.append(pn)
if len(pn_list)>0:
pnmean = np.mean(pn_list)
else:
pnmean=0
return pnmean
%matplotlib inline
df_tweets = pd.read_csv('tweets_nikkei.csv', names=['id', 'date', 'text', 'fav', 'RT'], index_col='date')
df_tweets = df_tweets.drop('text', axis=0)
df_tweets.index = pd.to_datetime(df_tweets.index)
df_tweets = df_tweets[['text']].sort_index(ascending=True)
means_list = []
for tweet in df_tweets['text']:
dl_old = get_diclist(tweet)
dl_new = add_pnvalue(dl_old, pn_dict)
pnmean = get_mean(dl_new)
means_list.append(pnmean)
df_tweets['pn'] = means_list
df_tweets = df_tweets.resample('D').mean()
df_tweets['pn'] = (df_tweets['pn'] - df_tweets['pn'].mean()) / df_tweets['pn'].std()
x = df_tweets.index
y = df_tweets.pn
plt.plot(x,y)
plt.grid(True)
df_tweets.to_csv('df_tweets_nikkei.csv')一応下のようなグラフは得られましたが、関係性は少なそうな印象です。
（ツイートの日付は３週間ずらしてプロットしてます。７日間で移動平均をかけてます。）
1.4. 東京の天気
気温、日照時間、雲量、降雨量、気圧、風速を取得。今回は結果的にあまり重要ではなかったのでグラフは割愛。
２．データセット
・説明変数：東京の感染者（７日後の増加人数、10万人当たり）
・説明変数：移動量、日経平均株価、日経電子版_マーケットのツイート、気温、日照時間、雲量、降雨量、気圧、風速（いずれも３週間ずらして７日間で移動平均処理）
３．解析
今回はランダムフォレストを用いて回帰分析しました。
コードは「化学のためのpythonによるデータ解析・機械学習入門」を引用しています。
（参照）https://datachemeng.com/data_analysis_in_chemistry_ohmsha/
from sklearn.ensemble import RandomForestRegressor
import matplotlib.pyplot as plt
import matplotlib.figure as figure
from sklearn import metrics
import numpy as np
import pandas as pd
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
#read_csv
df = pd.read_csv('dataset.csv',index_col=0)
#train_data
x_df = df.iloc[:, 1:]
y_df = df.iloc[:, 0]
x_train, x_test, y_train, y_test = train_test_split(x_df, y_df, test_size=0.25,random_state=1)
#auto_scaling
a_x_train = (x_train - x_train.mean()) / x_train.std()
a_y_train = (y_train - y_train.mean()) / y_train.std()
a_x_test = (x_test - x_train.mean()) / x_train.std()
#random_forest
rfr = RandomForestRegressor(n_estimators=100, random_state=0)
rfr.fit(a_x_train, a_y_train)
# テストデータで予測実行
estimated_y_test = rfr.predict(a_x_test)
# 特徴量の重要度を取得
feature = rfr.feature_importances_
# 特徴量の名前ラベルを取得
label = df.columns[1:]
# 特徴量の重要度順（降順）に並べて表示
indices = np.argsort(feature)[::-1]
for i in range(len(feature)):
print(str(i + 1) + " " +
str(label[indices[i]]) + " " + str(feature[indices[i]]))
# 特徴量の重要度の棒グラフ
plt.subplot(122, facecolor='white')
plt.title('Importance')
plt.bar(
range(
len(feature)),
feature[indices],
color='blue',
align='center')
plt.xticks(range(len(feature)), label[indices], rotation=90)
plt.xlim([-1, len(feature)])
plt.tight_layout()
# グラフの表示
plt.show()
# 説明変数の重要度
x_importances = pd.DataFrame(rfr.feature_importances_, index=x_df.columns, columns=['importance'])
x_importances.to_csv('rf_x_importances.csv')
##
# y を推定し、スケールをもとに戻す
estimated_y_train = rfr.predict(a_x_train) * y_train.std() + y_train.mean()
estimated_y_train = pd.DataFrame(estimated_y_train, index=y_train.index,columns=['estimated_y']) # Pandas の DataFrame 型に変換
R2 = r2_score(y_train, estimated_y_train)
print("R2: {}".format(R2))
# トレーニングデータの実測値 vs. 推定値のプロット
plt.rcParams['font.size'] = 18 # 横軸や縦軸の名前の文字などのフォントのサイズ
plt.figure(figsize=figure.figaspect(1)) # 図の形を正方形に
plt.scatter(y_train, estimated_y_train.iloc[:, 0], c='blue') # 実測値 vs. 推定値プロット
y_max = max(y_train.max(), estimated_y_train.iloc[:, 0].max()) # 実測値の最大値と、推定値の最大値の中で、より大きい値を取得
y_min = min(y_train.min(), estimated_y_train.iloc[:, 0].min()) # 実測値の最小値と、推定値の最小値の中で、より小さい値を取得
plt.plot([y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)],
[y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)], 'k-') # 取得した最小値-5%から最大値+5%まで、対角線を作成
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # y 軸の範囲の設定
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # x 軸の範囲の設定
plt.xlabel('actual y') # x 軸の名前
plt.ylabel('estimated y') # y 軸の名前
plt.show() # 以上の設定で描画
# トレーニングデータのr2, RMSE, MAE
print('r^2 for training data :', metrics.r2_score(y_train, estimated_y_train))
print('RMSE for training data :', metrics.mean_squared_error(y_train, estimated_y_train) ** 0.5)
print('MAE for training data :', metrics.mean_absolute_error(y_train, estimated_y_train))
# トレーニングデータの結果の保存
y_train_for_save = pd.DataFrame(y_train) # Series のため列名は別途変更
y_train_for_save.columns = ['actual_y']
y_error_train = y_train_for_save.iloc[:, 0] - estimated_y_train.iloc[:, 0]
y_error_train = pd.DataFrame(y_error_train) # Series のため列名は別途変更
y_error_train.columns = ['error_of_y(actual_y-estimated_y)']
results_train = pd.concat([estimated_y_train, y_train_for_save, y_error_train], axis=1)
results_train.to_csv('estimated_y_train.csv') # 推定値を csv ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください
# テストデータの推定
estimated_y_test = rfr.predict(a_x_test) * y_train.std() + y_train.mean() # y を推定し、スケールをもとに戻します
estimated_y_test = pd.DataFrame(estimated_y_test, index=x_test.index,
columns=['estimated_y']) # Pandas の DataFrame 型に変換。行の名前・列の名前も設定
# テストデータの実測値 vs. 推定値のプロット
plt.rcParams['font.size'] = 18 # 横軸や縦軸の名前の文字などのフォントのサイズ
plt.figure(figsize=figure.figaspect(1)) # 図の形を正方形に
plt.scatter(y_test, estimated_y_test.iloc[:, 0], c='blue') # 実測値 vs. 推定値プロット
y_max = max(y_test.max(), estimated_y_test.iloc[:, 0].max()) # 実測値の最大値と、推定値の最大値の中で、より大きい値を取得
y_min = min(y_test.min(), estimated_y_test.iloc[:, 0].min()) # 実測値の最小値と、推定値の最小値の中で、より小さい値を取得
plt.plot([y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)],
[y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)], 'k-') # 取得した最小値-5%から最大値+5%まで、対角線を作成
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # y 軸の範囲の設定
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # x 軸の範囲の設定
plt.xlabel('actual y') # x 軸の名前
plt.ylabel('estimated y') # y 軸の名前
plt.show() # 以上の設定で描画
# テストデータのr2, RMSE, MAE
print('r^2 for test data :', metrics.r2_score(y_test, estimated_y_test))
print('RMSE for test data :', metrics.mean_squared_error(y_test, estimated_y_test) ** 0.5)
print('MAE for test data :', metrics.mean_absolute_error(y_test, estimated_y_test))
# 誤差評価の保存
a_train = metrics.r2_score(y_train, estimated_y_train)
b_train = metrics.mean_squared_error(y_train, estimated_y_train) ** 0.5
c_train = metrics.mean_absolute_error(y_train, estimated_y_train)
a_test = metrics.r2_score(y_test, estimated_y_test)
b_test = metrics.mean_squared_error(y_test, estimated_y_test) ** 0.5
c_test = metrics.mean_absolute_error(y_test, estimated_y_test)
data1 = {"error": ["r^2","RMSE","MAE"], "train": [a_train, b_train, c_train], "test": [a_test, b_test, c_test]}
df1 = pd.DataFrame(data1)
df1.to_csv('gosa.csv')
# テストデータの結果の保存
y_test_for_save = pd.DataFrame(y_test) # Series のため列名は別途変更
y_test_for_save.columns = ['actual_y']
y_error_test = y_test_for_save.iloc[:, 0] - estimated_y_test.iloc[:, 0]
y_error_test = pd.DataFrame(y_error_test) # Series のため列名は別途変更
y_error_test.columns = ['error_of_y(actual_y-estimated_y)']
results_test = pd.concat([estimated_y_test, y_test_for_save, y_error_test], axis=1)
results_test.to_csv('estimated_y_test.csv') # 推定値を csv ファイルに保存移動量が大きく寄与する結果となりました、、
■学んだこと
・感染者、移動量などデータが公開されており、一般人も統計的なデータから考えることができる
・自粛が大事
この記事が気に入ったらサポートをしてみませんか？