一連のプロセスをひとまとめにしたPythonプログラムが以下のセルにあります.これは,同じフォルダーの中にあるdiagnostic_analytics.py
をそのままコピーしたものとなっています.これまでの説明と少し構成が異なる箇所もありますが(見つけてみて下さい),基本的に説明した内容で構成されています.
# -*- coding: utf-8 -*-
"""
The MIT License
Copyright 2021 WXBC,岐阜大学 吉野純
(C) 2021 WXBC,岐阜大学 吉野純
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
以下に定める条件に従い、本ソフトウェアおよび関連文書のファイル(以下「ソ
フトウェア」)の複製を取得するすべての人に対し、ソフトウェアを無制限に扱
うことを無償で許可します。これには、ソフトウェアの複製を使用、複写、変
更、結合、掲載、頒布、サブライセンス、および/または販売する権利、および
ソフトウェアを提供する相手に同じことを許可する権利も無制限に含まれます。
上記の著作権表示および本許諾表示を、ソフトウェアのすべての複製または重要
な部分に記載するものとします。
ソフトウェアは「現状のまま」で、明示であるか暗黙であるかを問わず、何らの
保証もなく提供されます。ここでいう保証とは、商品性、特定の目的への適合
性、および権利非侵害についての保証も含みますが、それに限定されるものでは
ありません。 作者または著作権者は、契約行為、不法行為、またはそれ以外で
あろうと、ソフトウェアに起因または関連し、あるいはソフトウェアの使用また
はその他の扱いによって生じる一切の請求、損害、その他の義務について何らの
責任も負わないものとします。
"""
# numpyをnpという別名でインポートします.
import numpy as np
# matplotlibをpltという別名でインポートします.
import matplotlib.pyplot as plt
# Seabornをインポートします.
import seaborn as sns
# datetimeは日時データを処理する際に便利なメソッドです.インポートします.
from datetime import date, datetime, timedelta
# pandasをpdという別名でインポートします.
import pandas as pd
# matplotlibで時系列図を作成するときには以下をインポートします
from pandas.plotting import register_matplotlib_converters
# これを登録しておきます.
register_matplotlib_converters()
# sklearn(scikit-learn)は機械学習関連のライブラリーです.インポートします.
from sklearn import linear_model
# 有意検定をするためにscipyのstatsというメソッドをインポートします.
import scipy.stats as stats
%precision 3
%matplotlib inline
#def 1.############################################
# 気象庁アメダスの気温の時系列データを読み込んで,
# データフレームに割り当てる関数
def readamedas(filename,skipline):
# pandasのread_csvというメソッドでcsvファイルを読み込みます.
# 引数として,
# [0]入力ファイル名
# [1]エンコーディング
# [2]読み飛ばす行数,
# [3]column名
# [4]datetime型で読み込むcolumn名
# [5]indexとするcolumn名
# を与える.
df=pd.read_csv(
filename,
encoding='Shift_JIS',
skiprows=skipline,
names=['date_time','temp','dummy1','dummy2'],
parse_dates={'datetime':['date_time']},
index_col='datetime'
)
return df
#def 2.############################################
# 東京電力の電力消費量の時系列データを読み込んで,
# データフレームに割り当てる関数
def readtepco(filename,skipline):
# pandasのread_csvというメソッドでcsvファイルを読み込みます.
# 引数として,
# [0]入力ファイル名
# [1]エンコーディング
# [2]読み飛ばす行数,
# [3]column名
# [4]datetime型で読み込むcolumn名
# [5]indexとするcolumn名
# を与える.
df=pd.read_csv(
filename,
encoding='Shift_JIS',
skiprows=skipline,
names=['date','time','power'],
parse_dates={'datetime':['date','time']},
index_col='datetime'
)
return df
#def 4.############################################
# 2つの時系列データからNaNを取り除いて,
# データフレームに割り当てる関数
def dataprocess0(df,name1,name2):
# 2つのデータのどちらかにNaNがあれば,その行を削除します.
df1=df.dropna(how='any')
# 処理した2つのname1とname2のデータを取り出し,
# 2つのデータ間の相関係数を計算します.
corr=df1[name1].corr(df1[name2])
# 処理したデータと相関係数を返します.
return df1,corr
#def 5.############################################
# 2つの時系列データから時系列図を作成する関数
def timeseries(df,name1,name2,filename):
# dfのインデックス(時間)をXとする
X=df.index
# dfのname1列を指定してデータを取り出し,numpy配列で値をY1に与える.
Y1=df.loc[:,[name1]].values
# dfのname1列を指定してデータを取り出し,numpy配列で値をY2に与える.
Y2=df.loc[:,[name2]].values
# 時系列図の大きさを指定
plt.figure(figsize=(20, 10))
# 1つ目(name1)のグラフを1行1列の1つ目に
ax1=plt.subplot(1,1,1)
# 2つ目(name2)のグラフのx軸を共有する
ax2=ax1.twinx()
# 1つ目(name1)の時系列
ax1.plot(X,Y1,color='blue',label=name1)
# 2つ目(name2)の時系列
ax2.plot(X,Y2,color='red',label=name2)
# グラフのタイトル
ax1.set_title("Timeseries:"+name1+" and "+name2)
# x軸のラベル
ax1.set_xlabel('Time')
# y軸(左側の第1軸)のラベル
ax1.set_ylabel('Temperature [degree C]')
# y軸(右側の第2軸)のラベル
ax2.set_ylabel('Electric Power [$10^6$kW]')
# 1つ目(name1)の凡例(左上に置く)
ax1.legend(loc='upper left')
# 2つ目(name1)の凡例(右上に置く)
ax2.legend(loc='upper right')
# 保存するファイル名
plt.savefig(filename)
# 図を閉じる
plt.close()
return
#def 6.############################################
# 2つの時系列データから散布図を作成する関数
def scatter(df,name1,name2,filename):
# scikit-learnの線形回帰モデルを呼び出す
clf = linear_model.LinearRegression()
# 説明変数Xにはname1を割り当てる(numpy配列)
X=df.loc[:,[name1]].values
# 説明変数Yにはname2を割り当てる(numpy配列)
Y=df.loc[:,[name2]].values
# Y=aX+bの予測モデルを作成する
clf.fit(X,Y)
# 回帰係数a
slope=clf.coef_
# 切片b
intercept=clf.intercept_
# 決定係数R2
r2=clf.score(X,Y)
# 文字列"Y=aX+b (R2=r2)"
equation=" y = "+str('{:.1f}'.format(slope[0][0]))+" x +"+str('{:.0f}'.format(intercept[0]))+" (R2="+str('{:.2f}'.format(r2))+")"
# 相関係数とその有意確率p-値を計算
corrcoef, pvalue = stats.pearsonr(np.ravel(X),np.ravel(Y))
# 散布図の大きさを指定
plt.figure(figsize=(8, 8))
# 散布図のプロット
plt.plot(X, Y, 'o')
# 散布図上に回帰直線を引く
plt.plot(X, clf.predict(X))
# 文字列"Y=aX+b (R2=r2)"を図の左上に置く
plt.text(np.nanmin(X), np.nanmax(Y), equation)
# グラフのタイトル
plt.title("Scatter diagram:"+name1+" and "+name2)
# x軸のラベル
plt.xlabel('Temperature [degree C]')
# y軸のラベル
plt.ylabel('Electric Power [$10^6$kW]')
# 保存するファイル名
plt.savefig(filename)
# 図を閉じる
plt.close()
return corrcoef, pvalue
#def 7.############################################
# 2つの時系列データから夏の期間(7/1から9/30まで)を取り出して,
# DataFrameに割り当てる関数
def dataprocess1(df,name1,name2):
# 夏の期間(6/1から9/30まで)のデータを取り出します.
df1=df[(df.index.date >= date(2017,7,1) ) & (df.index.date <= date(2017,9,30))]
# 処理した2つのname1とname2のデータを取り出し,
# 2つのデータ間の相関係数を計算します.
corr=df1[name1].corr(df1[name2])
# 処理したデータと相関係数を返します..
return df1,corr
#def 8.############################################
# 2つの時系列データから12~15時だけを取り出して,
# DataFrameに割り当てる関数
def dataprocess2(df,name1,name2):
# 14時のデータだけを取り出します.
df1=df[(df.index.hour >= 12)&((df.index.hour <= 15))]
# 処理した2つのname1とname2のデータを取り出し,
# 2つのデータ間の相関係数を計算します.
corr=df1[name1].corr(df1[name2])
# 処理したデータと相関係数を返します.
return df1,corr
#def 9.############################################
# 2つの時系列データから土曜日・日曜日のデータを取り除いて,
# DataFrameに割り当てる関数
def dataprocess3(df,name1,name2):
# 土曜日(weekday==5)と日曜日(weekday==6)を取り除く.
df1=df[(df.index.weekday != 5) & (df.index.weekday != 6) ]
# 処理した2つのname1とname2のデータを取り出し,
# 2つのデータ間の相関係数を計算します.
corr=df1[name1].corr(df1[name2])
# 処理したデータと相関係数を返します.
return df1,corr
#def 10.############################################
# 2つの時系列データからお盆や祝日のデータを取り除いて,
# DataFrameに割り当てる関数
def dataprocess4(df,name1,name2):
# お盆休暇(8/7から8/20まで)以外残す
df1=df[(df.index.date < date(2017,8,7) ) | (df.index.date > date(2017,8,20))]
# 海の日以外を残す
df1=df1[(df1.index.date != date(2017,7,17))]
# 山の日以外を残す
df1=df1[(df1.index.date != date(2017,8,11))]
# 敬老の日以外を残す
df1=df1[(df1.index.date != date(2017,9,18))]
# 秋分の日以外を残す
df1=df1[(df1.index.date != date(2017,9,23))]
# 処理した2つのname1とname2のデータを取り出し,
# 2つのデータ間の相関係数を計算します.
corr=df1[name1].corr(df1[name2])
# 処理したデータと相関係数を返します.
return df1,corr
###################################################
# ここから,メインプログラム
###################################################
#main 1.###########################################
# 気象庁AMeDAS(東京)の気温(℃)のcsvファイルの名前を指定します.
filename1='amedas2017.csv'
# csvファイルの上から5行目まではデータではないため呼び飛ばします.
skipline1=5
# 気象庁AMeDAS(東京)の気温(℃)のcsvファイルを読み込んで,
# pandasのデータフレーム(tepco)に割り当てる関数を呼び出します.
amedas=readamedas(filename1,skipline1)
# データフレーム(amedas)の中のdummy1とdummy2の列を削除する.
amedas=amedas.drop(['dummy1','dummy2'],axis=1)
# amedas
#main 2.##########################################
# 東京電力の電力消費量(10^6 kW)のcsvファイルの名前を指定します.
filename2='tepco2017.csv'
# csvファイルの上から3行目まではデータではないため呼び飛ばします.
skipline2=3
# 東京電力の電力消費量(10^6 kW)のcsvファイルを読み込んで,
# pandasのデータフレーム(tepco)に割り当てる関数を呼び出します.
tepco=readtepco(filename2,skipline2)
# tepco
#main 3.##########################################
# 2つのデータフレーム(amedasとtepco)を結合します.
data=pd.merge(amedas, tepco, how="outer", left_index=True, right_index=True)
# data
#main 4.##########################################
# データ前処理0: dataprocess0関数を呼び出して,
# どちらかにNaNがある行を取り除きます.
data0,corr=dataprocess0(data,'temp','power')
# data0
# 処理されたデータを用いて時系列図を作成します.
# ファイル名:timeseries0.png
timeseries(data0,'temp','power','timeseries0.png')
# 処理されたデータを用いて散布図を作成します.
# ファイル名:scatter0.png
corrcoef, pvalue=scatter(data0,'temp','power','scatter0.png')
# 処理された2つのデータの間の相関係数を表示します.
print("Dataprocess0... Corr=",corrcoef,"(p-value=",pvalue,")")
# csvファイルでエクスポート
data0.to_csv('data0.csv',sep=',')
#main 5.##########################################
# データ前処理1: dataprocess1関数を呼び出して,
# 夏の期間(7/1から9/30まで)だけを取り出します.
data1,corr=dataprocess1(data0,'temp','power')
# data1
# 処理されたデータを用いて時系列図を作成します.
# ファイル名:timeseries1.png
timeseries(data1,'temp','power','timeseries1.png')
# 処理されたデータを用いて散布図を作成します.
# ファイル名:scatter1.png
corrcoef, pvalue=scatter(data1,'temp','power','scatter1.png')
# 処理された2つのデータの間の相関係数を表示します.
print("Dataprocess1... Corr=",corrcoef,"(p-value=",pvalue,")")
# csvファイルでエクスポート
data1.to_csv('data1.csv',sep=',')
#main 6.##########################################
# データ前処理2: dataprocess2関数を呼び出して,
# 12~15時だけを取り出します.
data2,corr=dataprocess2(data1,'temp','power')
# data2
# 処理されたデータを用いて時系列図を作成します.
# ファイル名:timeseries2.png
timeseries(data2,'temp','power','timeseries2.png')
# 処理されたデータを用いて散布図を作成します.
# ファイル名:scatter2.png
corrcoef, pvalue=scatter(data2,'temp','power','scatter2.png')
# 処理された2つのデータの間の相関係数を表示します.
print("Dataprocess2... Corr=",corrcoef,"(p-value=",pvalue,")")
# csvファイルでエクスポート
data2.to_csv('data2.csv',sep=',')
#7.###############################################
# データ前処理3: dataprocess3関数を呼び出して,
# 土曜日・日曜日を取り除きます.
data3,corr=dataprocess3(data2,'temp','power')
# data3
# 処理されたデータを用いて時系列図を作成します.
# ファイル名:timeseries3.png
timeseries(data3,'temp','power','timeseries3.png')
# 処理されたデータを用いて散布図を作成します.
# ファイル名:scatter3.png
corrcoef, pvalue=scatter(data3,'temp','power','scatter3.png')
# 処理された2つのデータの間の相関係数を表示します.
print("Dataprocess3... Corr=",corrcoef,"(p-value=",pvalue,")")
# csvファイルでエクスポート
data3.to_csv('data3.csv',sep=',')
#7.###############################################
# データ前処理4: dataprocess4関数を呼び出して,
# お盆・祝日を取り除きます.
data4,corr=dataprocess4(data3,'temp','power')
# data
# 処理されたデータを用いて時系列図を作成します.
# ファイル名:timeseries4.png
timeseries(data4,'temp','power','timeseries4.png')
# 処理されたデータを用いて散布図を作成します.
# ファイル名:scatter4.png
corrcoef, pvalue=scatter(data4,'temp','power','scatter4.png')
# 処理された2つのデータの間の相関係数を表示します.
print("Dataprocess4... Corr=",corrcoef,"(p-value=",pvalue,")")
# csvファイルでエクスポート
data4.to_csv('data4.csv',sep=',')
課題1.冬季のデータを使って気温と電力消費量の関係性を調べてみよう.
# 上のプログラムをベースとして,解答プログラムをここに入力して下さい
課題2.2020年のデータを入手して,コロナ禍で関数dataprocess1
,関数dataprocess2
,関数dataprocess3
,関数dataprocess4
の仮説が正しいかどうかを調べてみよう.
# 上のプログラムをベースとして,解答プログラムをここに入力して下さい