東電消費電力データのベイズ曲線回帰
今回は東電の消費電力データにベイズ曲線回帰を適用してみました.
学習データは東電が公開しているcsvデータで,
入力は時刻 [時],出力は消費電力 [万kW]です.
(ほんとは気温とか,平日or休日かのダミー変数も入力しないと消費電力予測には使えないけど...)
理論式はPRML第1章を参考に以下のPythonコードを作成しました.
#coding:utf-8 ''' Purpose:Bayesian Curve_fitting ''' import numpy as np import input_csv import matplotlib.pyplot as plt # ////////////////////////// # ----- set data_set ----- # ///////////////////////// path = 'C:\Documents and Settings\User\My Documents\python_csv' read_file = path+'\\juyo-j.csv' head_num = 5 ; req_col = [2] # //////////////////////////// # ----- set parameters ----- # //////////////////////////// M = 10 alpha = 0.001 beta = 11.1 # ////////////////////////////// # ----- defined functions ----- # ////////////////////////////// def scale(X): # For Normalization mu = np.mean(X, axis=0) sigma = np.std(X, axis=0,ddof=1) X = (X-mu)/sigma return X,mu,sigma def phi(x,M): # This is one-column vector. return np.matrix([x**i for i in range(0,M+1)]).reshape(M+1,1) def S_inv(xlist,alpha,beta,M): tmp1 = sum([phi(x,M)*phi(x,M).T for x in xlist]) tmp2 = alpha*np.identity(M+1)+beta*tmp1 return np.matrix(tmp2) def mean(phi_x,S,xlist,tlist,beta,M): tmp_list = list(np.vstack((xlist,tlist)).transpose()) return beta*phi_x.T*S*sum([phi(xi,M)*ti for [xi,ti] in tmp_list]) def variance(phi_x,S,beta): return beta**(-1.)+phi_x.T*S*phi_x def main(): #>>----- get data_set ----- data = np.array(input_csv.method_01(read_file,head_num,req_col)) tlist,mu,sigma = scale(data[:,0]) xlist = np.array(range(0,len(tlist)))*0.1 #>>----- learn mean and variance ----- S = S_inv(xlist,alpha,beta,M).I xs = np.linspace(0, xlist[-1], 500) means = np.array([mean(phi(xi,M),S,xlist,tlist,beta,M)[0,0] for xi in xs]) s = np.sqrt([variance(phi(xi,M),S,beta)[0,0] for xi in xs]) plt.title('2011/4/5') plt.xlabel('Time[H]',style='italic',fontsize=15) plt.ylabel('power consumption [e4 kW]',style='italic',fontsize=15) plt.plot(xlist,tlist*sigma+mu,'ro') plt.plot(xs,means*sigma+mu,'b-') plt.plot(xs,(means+s)*sigma+mu,'g--') plt.plot(xs,(means-s)*sigma+mu,'g--') plt.xlim(0.0, xlist[-1]) plt.xticks(xlist, range(0,len(tlist))) plt.grid(True) plt.show() if __name__ == "__main__": main()
で,結果を描写したものがこれです.(赤点が学習データで緑の破線は分散による幅)
コード上にある'input_csv'モジュールは前回の記事で紹介したもので,
東電のcsv(juyo-j.csv)でヘッダが5行,用いた消費電力データが3列目だったので
"head_num = 5 ; req_col = [2]"としています.
なお,用いた多項式関数の次数Mやハイパーパラメータalpha,betaの値はテキトーです.
てか,東京電力電力供給状況APIで過去の消費電力取得できるけど,気温はどー取得すればいいんやろ?