東電消費電力データのベイズ曲線回帰

今回は東電の消費電力データにベイズ曲線回帰を適用してみました.

学習データは東電が公開している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で過去の消費電力取得できるけど,気温はどー取得すればいいんやろ?