プログラマーじゃない人が始める人工知能: ホールドアウト検証を使って最適な近似モデルを選択する
バグでベクトルに同じ値が入ってしまって、そのバグとりにめっちゃ時間がかかっていたのは内緒。
「プログラマーじゃない人が始める人工知能」の各記事はこちら。
【目次】
データを分割する
前回の記事で、架空会社のアクセス数はweek3.5以降くらいから急激に伸びていることを知ったオライリー先生と私。
直線の近似モデルではweek3.5以降のデータを上手く適合できていないのが見てとれるため、"いつの時点で100,000件を超えるのかを予測する"という目的を達成できそうにない。
なので、
- 次数を増やした近似モデルを複数用意する
- 各モデルの精度を比較して最適な近似モデルを選択する
という作業を実施していく。
まずはデータの準備だ。
最初にweek3.5以前のデータを捨ててしまおう。week3.5以降にサービスが急拡大しているようなので、それ以降の伸びを予測するためには、以前のデータは不要なのだ。
# 3.5週目からのデータを抽出 inflection = 3.5*7*24 xa = x[inflection:] ya = y[inflection:]
変数inflectionには、3.5週間を時間に直した値を入れる。(データの単位が時間なので) そして、SciPyの記法"[数値:]"を使って、数値より後の値のみを格納した配列を作成する。ちなみに、"[:数値]"と書けば数値より前の値のみを格納した配列を作成できる。
そして、ホールドアウトの準備にとりかかる。
既に手元にあるデータのうち、ある一定の割合をテストデータとし、残りのデータを訓練データとする。訓練データのみを用いてモデルを構築し、モデルで予測したデータをテストデータと比較してモデルの精度を検証する。こうすれば、手元にあるデータだけを使ってモデルの精度を検証できる。
1次関数の近似モデルの他に、2次、3次、10次、100次と次数を増やした近似モデルを追加して、その中から精度高いモデルを検証、選択する。
データのうち、30%をテストデータ。残り70%を訓練データとしよう。3割くらいをテストデータにするのが一般的のようだ。
# データをホールドアウトして配列を作成 frac = 0.3 split_idx = int(frac*len(xa)) shuffled = sp.random.permutation(list(range(len(xa)))) test = sorted(shuffled[:split_idx]) train = sorted(shuffled[split_idx:])
split_idxに3.5週以降のデータの値の数に30%を乗算した値をいれる。ちなみに、データの全数は735、3.5週以降の数は147、そのうち30%は44.1だ。int()で整数に直すと44。
shuffledには3.5週以降のデータをSciPyのランダム関数でシャッフルした配列を入れる。
そして、シャッフルした配列から44個の値をtestに、残りをtrainに入れる。
ホールドアウト検証する
データの準備が整ったので、各次数の近似モデルを構築しよう。polyfitにはtrainの配列のみをわたす。
# 各次数のモデル関数を作成 fbt1 = sp.poly1d(sp.polyfit(xa[train],ya[train],1)) fbt2 = sp.poly1d(sp.polyfit(xa[train],ya[train],2)) fbt3 = sp.poly1d(sp.polyfit(xa[train],ya[train],3)) fbt10 = sp.poly1d(sp.polyfit(xa[train],ya[train],10)) fbt100 = sp.poly1d(sp.polyfit(xa[train],ya[train],100))
polyfitでは最後の引数で次数を指定する。1なら1次、10なら10次。関数ごとpoly1dの引数にわたしてあげると、近似モデルの式をさっくり作ってくれる。
これらの近似モデルにアクセス件数を予測させ、実際のアクセス件数と比較することで、モデルの精度を検証したい。近似モデルの予測値と実際の値の誤差を定義しよう。
# 実際のデータyとモデル関数fにxを代入した際に算出される値の誤差 def error(f,a,b): return sp.sum((f(a)-b)**2)
f(a)がモデルが予測した値、bが実際の値だ。その差を計算して、二乗すれば、絶対値を求められる。
さあ、これで準備万端。
各予測モデルのアクセス件数と実際のアクセス件数を出力しよう。
for f in [fbt1,fbt2,fbt3,fbt10,fbt100]: print("Error d=%i: %f" %(f.order, error(f, xa[test],ya[test])))
モデル毎にprintしていけば初学者にもわかりやすいのだが、短く美しいコードを目指すオライリー先生がループ文を使いたいそうなので、その意向に従う。[fbt1...]と変数の分だけループしなさいよと宣言。次の行がループ内容。"Error..."文に続けて、fbt1～100の変数をerrorに渡して計算された値を表示する。結果はこんな感じだ。
Error d=1: 8458331.987260 Error d=2: 7855217.156308 Error d=3: 8007532.861243 Error d=10: 9746256.095559 Error d=53: 13058974.793724
当該データの量では100次の近似モデルが構築できず、次数が53になっているが、まあ、これはスルーしよう。あと、"RankWarning: Polyfit may be poorly conditioned"なる警告が出ているが、これも、まあ、スルーしよう。(※原因はイマイチ分からず。直し方をご存じの方がいらっしゃればコメントください)
近似モデルを用いて計算した値と実際の値の差が小さければ小さいほど、精度が高いといえる。結果を見ると、次数2のモデルの誤差が一番小さい。乱数が使われているので、実行の度に値が変化するが、大抵の場合、次数2が優秀だ。
最適な近似モデルを用いて予測する
というわけで、次数2の近似モデルが一番良さそうなことが分かった。
次数2の近似モデルを使用して、"いつの時点で100,000件を超えるのかを予測"しよう。100,000という値を代入して方程式を解けばよいのだが、Pythonではfsolveという関数を使う。
# fsolveをインポート from scipy.optimize import fsolve # 次数2のモデル関数を使って100,000アクセスの時点を予測 reached_max = fsolve(fbt2-100000,800)/(7*24) print("100,000 hits/hour expected at week %f" % reached_max[0])
プログラムを実行すると、大体、week9.8くらいというメッセージを得られる。
ついに結論を得た!
架空会社はおよそweek9.8くらいの時点で100,000件のアクセスに到達する。これで一安心だ。架空会社の情報システム部は胸をなでおろし、"week9.8まで"を目標にインフラの増強計画を進めていけばよい。
まあ、今がweek4が終わったところだから、時間はあと1か月ちょっとしかないんだけど、たぶん何とかなるよね。
今回書いたコードは以下の通り。Matplotlibを使って描画するところは除外。
# scipyとfsolveをインポート import scipy as sp from scipy.optimize import fsolve # テストデータを読み込み data = sp.genfromtxt("C:\Users\kazuma.tadaki\Desktop\Python_Sample_Programs\web_traffic.tsv", delimiter="\t") # 次元毎にデータを分割 x = data[:,0] y = data[:,1] # yの無効な値を除外 x = x[~sp.isnan(y)] y = y[~sp.isnan(y)] # 3.5週目からのデータを抽出 inflection = 3.5*7*24 xa = x[inflection:] ya = y[inflection:] # データをホールドアウトして配列を作成 frac = 0.3 # 30%を使ってテストする split_idx = int(frac*len(xa)) shuffled = sp.random.permutation(list(range(len(xa)))) test = sorted(shuffled[:split_idx]) train = sorted(shuffled[split_idx:]) # 各次数のモデル関数を作成 fbt1 = sp.poly1d(sp.polyfit(xa[train],ya[train],1)) fbt2 = sp.poly1d(sp.polyfit(xa[train],ya[train],2)) fbt3 = sp.poly1d(sp.polyfit(xa[train],ya[train],3)) fbt10 = sp.poly1d(sp.polyfit(xa[train],ya[train],10)) fbt100 = sp.poly1d(sp.polyfit(xa[train],ya[train],100)) # 実際のデータyとモデル関数fにxを代入した際に算出される値の誤差 def error(f,a,b): return sp.sum((f(a)-b)**2) for f in [fbt1,fbt2,fbt3,fbt10,fbt100]: print("Error d=%i: %f" %(f.order, error(f, xa[test],ya[test]))) # 次数2のモデル関数を使って100,000アクセスの時点を予測 reached_max = fsolve(fbt2-100000,800)/(7*24) print("100,000 hits/hour expected at week %f" % reached_max[0])