2020-04-25

ARIMA・SARIMAを用いてCOVID-19の感染者数を予測する

技術

はじめに

最近、新型コロナウイルスが世間を騒がせており、研究者が続々論文を発表したり、SIGNATEなどコンペサイトでも感染者数の予測コンペが行われるなど、様々な取り組みが行われています。
ここでは、時系列データ分析手法の一つであるARIMAとSARIMAを用いて日本の感染者数の予測を行いました。

ARIMAとSARIMAについて

ARIMAやSARIMAは複数の時系列データ分析手法を組み合わせたモデルになっています。
モデルの関係について、一般的に以下のような表で表されます。

f:id:resweater:20200424113436j:plain — （出典：未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge）

ARIMA

ARモデル、MAモデル、ARMAモデルはいずれも定常過程を対象とした時系列モデルです。定常過程に従う時系列データは、一定の周期で同程度の変動をしているデータのことです。しかし、常に一定のデータは滅多になく、非定常過程のデータであることが多いため、非定常過程の時系列データに対応した手法を用いる必要があります。
ARIMAモデルは非定常過程に対応したモデルの一つです。 ARIMAモデルは前後のデータ間の差分dを定義しており、ARモデルに関する次数p、MAモデルに関連する次数qと合わせてARIMA(p, d, q)と表すことができます。

SARIMA

SARIMAはARIMAモデルに長期的な季節変動を取り入れたモデルになります。実際の時系列データでは季節に合わせた周期性があるデータも多いため、それらを考慮したモデルとなっています。

予測で用いるデータ

COVID-19の感染者情報が載っているサイトはいくつかありますが、今回は以下のサイトからcsvファイルをインポートして用いました。

gis.jag-japan.com

読み込んだcsvファイルを加工し、1/15~4/20の日本全体の感染者数を抽出しました。

f:id:resweater:20200424115556p:plain — 1/15 ~ 4/20 1日あたりの全国感染者数

コレログラムと季節成分分解

元のデータと元のデータから時間をずらしたデータとの相関のことを自己相関といい、これを表したグラフをコレログラムと言います。コレログラムを見ることで、データの周期性を調べることができます。
グラフを見た感じ、ほとんど周期性はなく、時間が経つにつれて相関性が小さくなっていることが分かります。

f:id:resweater:20200424120326p:plain — コレログラム

次に、季節成分分解をします。これにより、トレンド・季節性・残差に分解することができ、データにどのようなトレンド性があるのかなどについて知ることができます。
グラフを見ると、時間の後半にかけてトレンドの数値が大きくなっていることが分かります。

f:id:resweater:20200424121033p:plain — 季節成分分解（上から順に元データ・トレンド・季節性・残差）

ARIMAでの予測

今回のデータは非定常過程のデータとみられることから、まずARIMA(p, d, q)を構築して予測してみることにしました（本当はデータが定常か非定常かについて検定をする必要がありますが、今回は省略しました、、、もっと勉強しなきゃ、、、）。
上記した通り、ARIMAにはp, d, qの3つのパラメータがあります。このうち、p, qはある程度自動推定してくれる機能があり、以下の通り実行します。

params = sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc')
print(params)

すると、以下のような結果が得られます。このうち、'aic_min_order': (4, 2)がp, qのパラメータになります。

{'aic':              0            1            2
0          NaN  1055.125479  1054.226782
1  1055.344019  1054.737810  1055.721728
2  1055.676846  1055.374114  1058.681375
3  1056.655811  1056.810369          NaN
4  1057.383730  1057.394580  1027.192160, 'aic_min_order': (4, 2)}

次に、差分dは自分で決めなければいけません。今回は差分回数は1で決め打ちします。差分データは以下の通り計算します。

diff = ts.diff()
diff = diff.dropna()

これでp, d, qの3つのパラメータが揃ったので、ARIMAモデルを構築します。tsは時系列データです。

from statsmodels.tsa.arima_model import ARIMA
arima_model = ARIMA(ts, order=(4,1,2)).fit(dist=False)

ここで、モデルから推測される値と実測値の差（残差）を表示してみます。後半になるにつれ残差が大きいことが分かります。

f:id:resweater:20200425180558p:plain — ARIMAモデルの残差

次に、構築したARIMAモデルを用いて予測をしてみます。青い線が1/15 ~ 4/20の実際の値、オレンジ色の線が4/20 ~ 5/10の予測値になります。かなり大きくぶれていて、あまり上手く予測できていなさそうなのが分かります。

f:id:resweater:20200425180905p:plain — 4/20 ~ 5/10のARIMAモデルによる予測結果

SARIMAでの予測

SARIMA(p, d, q)(P, D, Q)[s]はARIMAのモデルのパラメータに加え、季節調整のsというパラメータとP, D, Qというパラメータが必要になります。
ここで、今回のデータはコレログラムでは周期性が見られませんでしたが、ドメイン知識として毎週月曜日は検査数が減るため感染者数も減るという事が分かっています。そこで、今回は7日ごと周期性があると仮定して、s=7でモデル構築を行いました。
P, D, Qは本来はグリッドサーチなどで最適な値を求める必要がありますが、今回はいずれも1で決め打ちしました。

sarima_model = sm.tsa.SARIMAX(ts, order=(4,1,2), seasonal_order=(1,1,1,7)).fit()

残差をプロットしてみると以下のようになりました。ARIMAモデルより振り幅が小さいことが分かります。

f:id:resweater:20200425181553p:plain — SARIMAモデルの残差

構築したSARIMAモデルを用いた予測結果は以下の通りです。ARIMAモデルと同じく、青い線が1/15 ~ 4/20の実際の値、オレンジ色の線が4/20 ~ 5/10の予測値になります。ARIMAモデルよりもそれっぽい予測が出来ていることが分かります。

f:id:resweater:20200425181825p:plain — 4/20 ~ 5/10のSARIMAモデルによる予測結果

まとめ

時系列データに対してARIMAモデルとSARIMAモデルを用いて予測を行うことで、それっぽい結果を得ることができました。ただ、この結果はあくまで時系列データしか使っていないため、実際は様々な外的要因によって変化していくと思われます。
今回、初めて時系列データのモデル構築に挑戦したので、おかしいところがあればぜひご指摘下さい。あと、グラフの時間軸が上手くプロットできなかったのが今後の課題です。。。

プログラム

分析で用いたプログラムはGitHubで公開しています。 github.com

参考

2019-12-30

2019年を振り返る

今年も様々な取り組みをしてきたので、振り返りつつまとめたいと思います。

研究

今年は2つの研究に取り組みました。
1つ目の研究はEEGから画像を生成した上でCNNで学習を行うという研究で、9月には国際学会 (SOMET2019) で発表を行いました。

resweater.hatenablog.com

2つ目の研究はCNNのハイパーパラメータ最適化に関する研究で、現在ジャーナルに投稿中です。
また、この研究の過程で得た知識を記事にしました。

resweater.hatenablog.com

qiita.com

データ分析・機械学習に関するあれこれ

コンペ

今年はいくつかコンペに参加し、中でもSIGNATEで行われたAI創薬コンペでは銀圏に入ることができました。

resweater.hatenablog.com

データ分析・機械学習を教える

現在、知能システムに関する研究室に所属しています。
TAとして後輩にデータ分析や機械学習の基礎を教える機会が与えられたため、試行錯誤しながら指導しました。

resweater.hatenablog.com

好きなアーティストの分析

気分転換をしたいときに行った分析として、TrySailというグループの分析を行っていました。
Spotifyに登録されている音楽情報を基に統計分析を行ったり、TrySailの夏川椎名さんのブログを学習して文章生成を行ったりしました。

resweater.hatenablog.com

github.com

その他

就活

大学院を20卒として修了する予定のため、1~3月を中心に就活を行いました。
結果として、データ分析事業を行っている企業から内定をいただくことができました。

resweater.hatenablog.com

宇宙の話

岩手で毎年行われているいわて銀河フェスタに参加し、世界最先端の宇宙に関する研究に触れることができました。

resweater.hatenablog.com

LT会

大学内で定期的に行われているLT会に参加し、発表を行いました。

resweater.hatenablog.com

まとめ

今年も様々な取り組みを行うことができました。
特に研究で2本論文を書けたこと、データ分析事業を行っている企業から内定をいただけたことは大きな成果でした。
来年も多くの取り組みを行ってブログでもアウトプットしていきたいと思います。

2019-11-21

SIGNATE「AI創薬：薬物動態予測」コンペで17位でした

技術

はじめに
創薬コンペ概要
試した手法
- メインモデル
- その他のあれこれ
結果
おわりに

はじめに

SIGNATEで行われた「AI創薬：薬物動態予測」コンペ（以下，創薬コンペ）に参加し，17位で上位3%に入る結果となりました．

signate.jp

本記事では，コンペで行った手法やその振り返りをまとめます．
ただし，コンペには参加規約*1があります．今回，自分は入賞者ではないので，非入賞者としての参加規約に反しない範囲でまとめたいと思います．

創薬コンペ概要

本コンペでは，新薬開発で重要な因子の一つである化合物の薬物動態パラメータを予測しました．
用いるデータとしては，様々な化合物の情報とそのスクリーニング試験における実験データが与えられました．ただ，どの特徴量も説明は非公開となっており，ドメインを生かした予測は難しくなっていました．

評価方法は決定係数「R2」でした．
R2は精度が高いほど大きな値となる評価指標です．
f:id:resweater:20191020140223p:plain

試した手法

メインモデル

f:id:resweater:20191021140058p:plain
コンペでは主に図に示したような流れで得た結果をサブミットました．
特にData Cleaningはするかしないかで大きな差が出ました．これについては以下のページを参考にさせてもらいながら処理を行いました．
naotaka1128.hatenadiary.jp

モデルの構築に当たっては，精度が出そうなモデルを中心に片っ端から手法を試す方法を取りました．
様々なモデルを試す中で，結果の傾向を見ながら複数のモデルの結果でStackingを行い，その結果をサブミットしました．

その他のあれこれ

上記以外の方法として，以下のような方法を試してみました．

データをカテゴリカルデータと数値データで分け，それぞれにあった手法で学習して統合
→ あまり精度が出なかったのと，並行してチームメンバーがやっていた方法が上手くいったためやめました．もう少し工夫してやってみたかった方法でした．

寄与度を基にデータを分けて学習して統合
→ 寄与度に一定の傾向が見られたため，その辺りを工夫すれば上手くいくかなと思い試しましたが，あまり良い結果は得られませんでした．

結果

f:id:resweater:20191020143840p:plain

Public LBからPrivate LBで若干精度は上がりましたが，順位はほぼ変わらず17位という結果になりました．

おわりに

最終的な参加チームは729チームだったので，かなり上位の方に食い込むことができました．ただ，コンペ終了数日前は9位まで上がり，目標は10位以内だったので惜しい結果となりました．
今回得た前処理などの手法を生かしてまた別のコンペにも参加していきたいと思います．

*1:https://signate.jp/competitions/168#terms

2019-09-27

初めて国際学会に参加してきた話

技術イベント

先日，初めて国際学会に参加してきたのでそれについてまとめます．

学会について

参加した学会

今回はマレーシアのクチンで行われた SOMET2019 (THE 18TH INTERNATIONAL CONFERENCE ON INTELLIGENT SOFTWARE METHODOLOGIES, TOOLS, AND TECHNIQUES) に参加・発表してきました．
SOMETは名前の通り知能システムやそれに関連するソフトウェア，システムなどを対象とした学会です．例年は3日間の日程で行われますが，今年は2日間に短縮され，3日目はクルージングツアーが組まれていました．

f:id:resweater:20190927153110j:plain

SOMET全般を通して

自分は学会というものに対する参加自体が初めてだったため，一つ一つがとても新鮮に感じました．
特にcoffee breakは様々な教授や学生と話すことができ，楽しく過ごすことができました．
また，keynote speakerのAIと5Gの今後に関する話や優秀な方々の発表は聞いてるだけで刺激を受けました．

自分の発表

自分はCNNを用いた脳波 (EEG) からの感情認識に関する研究を行っています．今回の学会では，"Emotion Recognition by Convolutional Neural Network Based on EEG-Images Plotting Time Series Data" というタイトルで発表してきました．
これは，EEGを一度時系列データとしてプロットして画像化し，その画像に対してCNNを用いて感情認識を行うというものです．

自分の発表は2日目の最終セッションでしたが，1日目に話した優秀な方から聞きに行くねと言われたりしてとても緊張しました……
また，学会が2日間に短縮された関係で発表時間も短縮させられたため，発表直前まで調整をし続けることになり大変でした．
実際の発表はすんなり話すことができました．しかし，質疑応答で難しい質問が飛んできた際，上手く英語で答えることができず，自分の力不足を感じました．

f:id:resweater:20190927153201j:plain

クチンという町について

学会が行われたクチンはリゾート地区であったため，観光していて楽しかったです．特に毎晩川沿いで行われるウォーターショーは見ごたえがあり，町の活気を感じました．
また，3日目に行われたクルーズでは，ちょっと町の中心部から離れると見えてくるクチンに暮らす人達の生活が垣間見えて興味深かったです．

f:id:resweater:20190927153318j:plain

まとめ

初めての学会でしたが多くの体験をし，刺激を受けることができました．
特に学会の質疑応答をはじめとして，いくつかの場面で自分の英語が上手く伝わらないことがあったので，研究を継続しつつ英語力も鍛えて，またどこかの学会でリベンジしたいと思いました．

2019-08-28

データ分析・機械学習について学ぶゼミを運営した話

技術

私は現在，知能システムに関する研究室に所属しています．担当教授と話す中で，後輩（3年生）のデータ分析に関して勉強する演習ゼミの運営を院生メインで行うことになりました．
この記事では，現状これといって正解のないデータ分析や機械学習の教え方について，どうやって教えていったのか，ゼミを運営していく中で感じたことなどについてまとめます．

現状と目標
- 現状
- 目標
実際の取り組み
やってみて見えた課題と良かったところ
終わりに

現状と目標

現状

ゼミを運営していく上で，まず教える相手の現在の理解レベルについて把握する必要があります．弊学の場合，一般的な3年生では以下のような学習レベルになっています．

C，Javaあたりは必修である程度学習済みだが，Pythonは学んでいない
統計学や線形代数，離散数学の講義は受けている
データ分析・機械学習は完全に未経験

データ分析について勉強していくことを考えると少し厳しい状況ですが，世の中のデータ分析を学んでみたい人というのは割とこのレベルが多い気がします．．．

目標

演習ゼミは前期に週に二回，合計30回程度あるため，この中でできることとして以下のような目標を立てました．

一般的な機械学習手法や関連技術について，名前とどのように使うかについてざっくり理解する
数式レベルまでの理解は求めない（一応説明は行う）
今後卒業研究を行う上で必要となる考え方について理解する

実際の取り組み

各回のゼミでは，基本的に以下の手順で運営を行いました．

教員や院生から技術について説明
演習
次回までの課題提示

教えた内容については，以下の順番で30回程度に分けて教えました．機械学習手法と関連技術を同時に紹介していくという手順です．
演習や課題として取り組んでもらったデータ分析では，UCI Machine learning repositoryやSIGNATEの練習問題を用いて分析を行ってもらいました．

Jupyter notebookの使い方・Pythonになれる・線形回帰
Support Vector Machine (SVM)・データの見方
ロジスティック回帰・Cross Validation
決定木・ランダムフォレスト
単純パーセプトロン・前処理
Neural Network (NN) フルスクラッチチャレンジ
Convolutional Neural Network (CNN)
教師なし学習について
興味のあるテーマ調査・院生の研究紹介・外部講師講演など

やってみて見えた課題と良かったところ

課題

1講義あたりの時間が短くてなかなか演習までいけず，時間内に収めるのが難しい
教えるための資料作成に時間が取られ，自分の研究が止まりがちになる
時系列データやRNN，LSTM，半教師あり学習，強化学習については教えられていない
一度単純パーセプトロンやNNの説明をしていても，NNのフルスクラッチを講義内で行うことは難しい

良かったところ

一般的な手法や関連技術について知ってもらい，一回は使ったことがあるレベルまで上げることができた
データ分析に関する知識をある程度持ってもらった上で，今後取り組みたい研究テーマについて考えてもらうことができた
毎回資料を作成するため，自分たちの復習になり，新たな発見なども得られる
機械学習を人に教えるということを学ぶことができる

終わりに

半年間ゼミを運営してみて，とりあえず目標は達成することができたと考えています．
とにかくデータ分析が楽しいと感じてもらえるよう，SIGNATEにあるような現実的なデータを取り入れたりして演習を行いました．
また，演習を取り入れたことで，個人個人の理解度を把握することができ，受け身にならず考えながら取り組んでもらえるようにしました．
ゼミを受けた3年生からゼミをやってみてのアンケートも行いましたが，全体的に好印象のようで安心しました……
後期以降もアンケートや夏休み中に各自行ったことなどを参考にしつつ，試行錯誤しながら運営していこうと考えています．

2019-08-26

いわて銀河フェスタでEHTや天体観測におけるデータの前々処理について聞いてきました

いわて銀河フェスタについて

8月24日 (土) ，国立天文台水沢キャンパスでいわて銀河フェスタ2019いわて銀河フェスタ2019が行われました．このイベントでは国立天文台の特別開放が行われ，今年話題になったブラックホールに関する研究 (EHT: Event Horizon Telescope) に携わった研究者方の講演会や座談会なども行われました．
イベントでは，天文学に携わる世界屈指の研究者の方々のお話を伺ったり交流したり~~できる一方，場所が水沢という田舎の地な上に広報が上手くされていないため，悠々と気が済むまでお話をさせていただくことが~~できました．
この記事では，そこで見たものや聞いてきた話についてまとめます．

見たもの・聞いたもの

ILCについて

会場に着いてまず，国際リニアコライダー (ILC: International Linear Collider) *1の広報ブースを訪ねました．ILCは素粒子物理学の研究で必要となる次世代の直線型衝突加速器で，岩手～宮城に跨る北上山地の地下へ建設することを目指して誘致活動が行われています．
現在は学術会議や研究者，政府で様々な意見が出ており*2 *3，建設計画が進むか予断が許さない状況です．
伺った話によると，2020年にマスタープランが策定されるまで主に3つの意思決定の場があり，それぞれの場でどのような判断がされるかについて注視しているとのことでした．
自分が所属している大学でもILCの建設を見越したプロジェクトや研究の話はいくつも聞いたことがあり，今後どのように進んでいくのか注目しています．

スーパーコンピュータ「アテルイⅡ」

国立天文台水沢キャンパスには，天文学研究用スーパーコンピュータの「アテルイⅡ」が置かれています．昨年までは「アテルイ」が稼働していましたが，今は一新され「アテルイⅡ」が稼働しています．
今回のイベントでも公開されていたのですが，前述したILCのブースで話をしている間にスパコン室に入るための整理券が無くなってしまい，入ることができませんでした……
ただ，整理券が取れなかった人向けに外からチラ見せをさせてもらいました．

f:id:resweater:20190825231030j:plain — 天文学専用スーパーコンピュータ「アテルイⅡ」

相関器室でデータの前々処理について学ぶ

相関器室では，VERAプロジェクトにおける観測データの処理についてお話を伺うことができました．

VERAについて

VERAはVLBIと呼ばれる手法によって天体の距離や運動を高精度で記録するプロジェクトです*4．高精度に観測を行うため，水沢（岩手）・入来（鹿児島）・石垣島・小笠原の4ヶ所の観測局を組み合わせて観測を行うことで，直径2,300kmの望遠鏡と同程度の性能を発揮しています．
近年は韓国や東アジアにある観測局を組み合わせることで，さらに精度を上げようとしたりもしているそうです．

f:id:resweater:20190825233001j:plain — 観測で使われる20mアンテナ．手前の線は北緯39°8'を示している．

取り込むデータと利用までのプロセス

20mアンテナを使うことで，様々な電波を受信することができます．しかし，raw dataをそのまま使うことは難しく，大抵の場合は他の観測局で取得できたデータと合わせて処理を行った上で，ユーザや研究者にデータが渡されます．これら相関処理と呼ばれる前々処理を行っているのが相関器室だそうです．
先日話題になったブラックホールの発見に関する研究でも，取得したraw dataに対して相関処理を行った上で，データをスパースモデリングによってイメージング処理を行っているとのことでした．

f:id:resweater:20190826002552j:plain — 取得したデータが研究者に届くまで

相関器室での質問と回答

相関器室では研究員の方が案内をしてくださり，また建物の奥の方の部屋のためか人も少なかったため，多くの質問をさせていただくことができました．

Q: 取得したraw dataをそのまま使うことはないのか
A: 相関処理は共通して行う前々処理なので，基本的にはない．ただ，プロジェクトによってはそのままのデータを渡すこともある．

Q: ブラックホールの研究も同じことをしているのか．
A: そう．ブラックホールの研究も相関処理を行ったうえでイメージング作業を行っている．

Q: 最近は人工衛星が天体観測に悪影響を与えているというニュースもあるが*5，実際どのような影響があるのか．
A: 今行っているプロジェクトでは，取得する電波と人工衛星が発している電波の周波数帯が異なるため影響はあまりない．ただ，プロジェクトによっては周波数帯が重なってしまう場合もあり，一回でも人工衛星の電波を取得してしまうと電波強度が強すぎるためノイズとして処理することは難しく，データを捨てるしかなくなる．

Q: 機械学習やAIと呼ばれるような技術は使われているか．もしくは使われる余地はあるのか．
A: 今のところ機械学習は取り入れていない．しかし，運用における異常検知や，一定間隔で観測している天体の移動における異常差をDeep Learningなどで自動検知できれば，手作業で行っている部分が自動化できて良いと思う．

EHT研究者の講演と座談会

ブラックホール発見の研究プロジェクトに携わった研究者の方々の講演会と座談会も行われ，プロジェクト中の裏話などを聞くことができました．
小さい子供や高齢者の方も多かったので，質疑応答の時も緩い質問が多いのかなと思っていましたが，「2年前にアメリカで発表されたブラックホールの重力波発見と比べて今回の研究成果はどのような違いや意義があるのか」というような非常にレベルの高い質問も出てきて興味深かったです．
また，子供からの「ブラックホールに行ってみたいか」という質問に対して，本間先生が「寿命が分かっててもうすぐ死ぬというタイミングだったら，ブラックホール葬じゃないけど行ってみたい笑」と答えていたことが面白かったです．

終わりに

天文学は自分の専攻ではありませんが，間接的に関連している話や技術も多く，非常に興味深かったです．また，非常にレベルの高い研究者の方々とこれほど気楽に話すことのできるイベントはないと思うので，貴重な経験をすることがきました．

*1:ILC Project

*2:学術会議、次世代加速器「誘致支持せず」巨額建設費などネック

*3:【速報】ＩＬＣ、議論継続へ　政府、誘致表明至らず

*4:VERA

*5:SpaceXの人工衛星群、天体観測に悪影響国立天文台など懸念示す

2019-08-22

ハイパーパラメータ最適化ツールのHyperactiveを試す

技術

この記事では，ハイパーパラメータ最適化ツールのHyperactiveについてまとめています．

ハイパーパラメータ最適化について
Hyperactive について
Random ForestのハイパーパラメータをPSOで最適化
- 使うデータセット
- 最適化なしと最適化を行った上での分類結果
CNNのハイパーパラメータをRandom Searchで最適化
- 使うデータセット
- 最適化なしと最適化を行った上での分類結果
終わりに

ハイパーパラメータ最適化について

機械学習では人が手作業で調整する必要があるハイパーパラメータが存在します．ハイパーパラメータを調整することは高い精度を求める上で重要となりますが，手作業では限界があるため，多くの自動化手法が提案されています．
提案されていたり利用されている自動化手法やツールとして，以下のようなものがあります．

Grid Search
Optuna (TPE)
Meta-Heuristics Algorithm (e.g. PSO, ES)

Hyperactiveは，これまであまり提案されてこなかったMeta-Heuristics Algorithmやその他の手法を用いて最適化を行うことができます．

Hyperactive について

概要

Hyperactiveは機械学習や深層学習のハイパーパラメータを自動で最適化するAPIです．このAPIを使うことで，様々な機械学習手法のハイパーパラメータに対して進化戦略 (ES) や粒子群最適化 (PSO) などの手法を用いて最適化を行うことができます．
この記事を書いている時点 (2019/08/21) でver. 0.4.1.5であり，現在も開発が行われています．

github.com

使うことの出来る機械学習パッケージと最適化手法

Hyperactiveでは，以下の機械学習パッケージに対応しています．

scikit-learn
LightGBM
CatBoost
Keras

また，Hyperactiveは以下の通り多くの最適化手法を使うことができます．

Hill Climbing
Stochastic Hill Climbing
Tabu Search
Random Search
Random Restart Hill Climbing
Random Annealing
Simulated Annealing
Stochastic Tunneling
Parallel Tempering
Particle Swarm Optimization
Evolution Strategy
Bayesian Optimization

インストール

READMEにある通り，次のようにインストールします．

pip install hyperactive

Random ForestのハイパーパラメータをPSOで最適化

今回は例としてRandom Forestのハイパーパラメータ最適化を行いました．

使うデータセット

今回はscikit-learnにあるIrisデータセットを用いました．Irisは4次元の特徴量，3種類のラベルを持つアヤメのデータセットです．

最適化なしと最適化を行った上での分類結果

まずRandom Forestを用いてハイパーパラメータ最適化を行わずに結果を求め，その後HyperactiveのPSOを用いて最適化を行った上で結果を求めました．

最適化なし

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
import time

iris_data = load_iris()
X = iris_data.data
y = iris_data.target

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

t1 = time.time()
forest = RandomForestClassifier()
forest.fit(X_train, y_train)
t2 = time.time()

print('Train score: {}'.format(forest.score(X_train, y_train)))
print('Test score: {}'.format(forest.score(X_test, y_test)))
print("time: {}".format(t2-t1))

結果

Train score: 0.9809523809523809
Test score: 0.9111111111111111
time: 0.009917736053466797

PSOで最適化

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from hyperactive import SimulatedAnnealingOptimizer
from hyperactive import ParticleSwarmOptimizer
import time

iris_data = load_iris()
X = iris_data.data
y = iris_data.target

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

# this defines the model and hyperparameter search space
search_config = {
    "sklearn.ensemble.RandomForestClassifier": {
        "n_estimators": range(10, 100, 10),
        "max_depth": [3, 4, 5, 6],
        "criterion": ["gini", "entropy"],
        "min_samples_split": range(2, 21),
        "min_samples_leaf": range(2, 21),
    }
}

# start point
start_point = {
    "sklearn.ensemble.RandomForestClassifier.0": {
        "n_estimators": [30],
        "max_depth": [6],
        "criterion": ["entropy"],
        "min_samples_split": [12],
        "min_samples_leaf": [16],
    }
}

#Optimizer = SimulatedAnnealingOptimizer(search_config, n_iter=100, n_jobs=4)
Optimizer = ParticleSwarmOptimizer(search_config, 
                                    n_iter=10,  # number of iterations to perform
                                    metric="accuracy", 
                                    n_jobs=1, 
                                    cv=3, 
                                    verbosity=1, 
                                    random_state=None, 
                                    warm_start=start_point,  # Hyperparameter configuration to start from
                                    memory=True,  # Stores explored evaluations in a dictionary to save computing time
                                    scatter_init=False,  # Chooses better initial position by training on multiple random positions with smaller training dataset 
                                    n_part=10,  # number of particles
                                    w=0.5,  # interia factor
                                    c_k=0.5,  # cognitive factor
                                    c_s=0.9)  # social factor

# search best hyperparameter for given data
t1 = time.time()
Optimizer.fit(X_train, y_train)
t2 = time.time()
print("time: {}".format(t2-t1))

# predict from test data
prediction = Optimizer.predict(X_test)

# calculate accuracy score
score = Optimizer.score(X_test, y_test)

print("test accracy: {}".format(score))

結果と最適化されたパラメータ

start_point = {'sklearn.ensemble.RandomForestClassifier.0': {'n_estimators': [70], 'max_depth': [4], 'criterion': ['gini'], 'min_samples_split': [16], 'min_samples_leaf': [13]}}
train score: 0.9721780604133546
test score: 0.9555555555555556
time: 18.719361782073975

Hyperactiveにも各手法ごとハイパーパラメータが存在しますが，今回は初期値のままで行いました．最適化されたパラメータはstart_pointに格納されるようです．

CNNのハイパーパラメータをRandom Searchで最適化

次に，深層学習の例としてCNNのハイパーパラメータ最適化を行いました．今回はGPUを積んでいないノートPCで実験を行っており，CNNはそれ自体の計算時間が長いため，今回は最適化手法の中で比較的計算コストの低いRandom Searchを用いています．また，いずれもepoch=5，batch size=500とし，損失関数にはcategorical clossentropyを用いています．

使うデータセット

ここではKerasにあるMNISTデータセットを用いました．MNISTは784 (28×28) 次元の特徴量，10種類のラベルを持つ手書き数字の画像データセットです．

最適化なしと最適化を行った上での分類結果

最適化なし

import time
import numpy as np
from keras.datasets import mnist
from keras.utils import to_categorical
from keras.layers import Input, Dense, Activation, BatchNormalization, Flatten, Conv2D, MaxPool2D, Dropout
from keras.models import Model
from keras.models import Sequential

(X_train, y_train), (X_test, y_test) = mnist.load_data()

X_train = X_train.reshape(60000, 28, 28, 1)
X_test = X_test.reshape(10000, 28, 28, 1)

y_train = to_categorical(y_train)
y_test = to_categorical(y_test)


model = Sequential()
model.add(Conv2D(filters=32, kernel_size=(3,3), 
                activation='relu', input_shape=(28,28,1)))
model.add(MaxPool2D(pool_size=(2,2)))
model.add(Conv2D(filters=(32), kernel_size=(3,3), 
                activation='relu'))
model.add(MaxPool2D(pool_size=(2,2)))
model.add(Flatten())
model.add(Dense(30, activation='relu'))
model.add(Dropout(0.4))
model.add(Dense(10, activation='softmax'))

model.compile(optimizer='sgd',
                loss='categorical_crossentropy',
                metrics=['accuracy'])

t1 = time.time()
history = model.fit(X_train, y_train, 
                    epochs=5, 
                    batch_size=500, 
                    verbose=1,
                    validation_data=(X_test, y_test))
t2 = time.time()
print("time: {}".format(t2-t1))

score = model.evaluate(X_test, y_test, batch_size=32, verbose=0)
print('validation loss:{0[0]}\nvalidation accuracy:{0[1]}'.format(score))

結果

train score: 0.2430
test score: 0.3045
time: 146.26237177848816

Random Searchで最適化を行った結果

import time
import numpy as np
from keras.datasets import mnist
from keras.utils import to_categorical

from hyperactive import RandomSearchOptimizer, ParticleSwarmOptimizer

(X_train, y_train), (X_test, y_test) = mnist.load_data()

X_train = X_train.reshape(60000, 28, 28, 1)
X_test = X_test.reshape(10000, 28, 28, 1)

y_train = to_categorical(y_train)
y_test = to_categorical(y_test)


# this defines the structure of the model and print("time: {}".format(t2-t1))the search space in each layer
search_config = {
    "keras.compile.0": {"loss": ["categorical_crossentropy"], "optimizer": ["SGD"]},
    "keras.fit.0": {"epochs": [5], "batch_size": [500], "verbose": [2]},
    "keras.layers.Conv2D.1": {
        "filters": [32, 64],
        "kernel_size": [(3, 3)],
        "activation": ["relu"],
        "input_shape": [(28, 28, 1)],
    },
    "keras.layers.MaxPooling2D.2": {"pool_size": [(2, 2)]},
    "keras.layers.Conv2D.3": {
        "filters": [16, 32],
        "kernel_size": [(3, 3)],
        "activation": ["relu"],
    },
    "keras.layers.MaxPooling2D.4": {"pool_size": [(2, 2)]},
    "keras.layers.Flatten.5": {},
    "keras.layers.Dense.6": {"units": [30], "activation": ["relu"]},
    "keras.layers.Dropout.7": {"rate": list(np.arange(0.4, 0.8, 0.2))},
    "keras.layers.Dense.8": {"units": [10], "activation": ["softmax"]},
}

Optimizer = RandomSearchOptimizer(search_config, n_iter=10)

# search best hyperparameter for given data
t1 = time.time()
Optimizer.fit(X_train, y_train)
t2 = time.time()

print("time: {}".format(t2-t1))

# predict from test data
prediction = Optimizer.predict(X_test)

# calculate accuracy score
score = Optimizer.score(X_test, y_test)

print("test accracy: {}".format(score))

結果と最適化されたパラメータ

start_point = {'keras.compile.0': {'loss': ['categorical_crossentropy'], 'optimizer': ['SGD']}, 'keras.fit.0': {'epochs': [5], 'batch_size': [500], 'verbose': [2]}, 'keras.layers.Conv2D.1': {'filters': [32], 'kernel_size': [(3, 3)], 'activation': ['relu'], 'input_shape': [(28, 28, 1)]}, 'keras.layers.MaxPooling2D.2': {'pool_size': [(2, 2)]}, 'keras.layers.Conv2D.3': {'filters': [16], 'kernel_size': [(3, 3)], 'activation': ['relu']}, 'keras.layers.MaxPooling2D.4': {'pool_size': [(2, 2)]}, 'keras.layers.Flatten.5': {}, 'keras.layers.Dense.6': {'units': [30], 'activation': ['relu']}, 'keras.layers.Dropout.7': {'rate': [0.4]}, 'keras.layers.Dense.8': {'units': [10], 'activation': ['softmax']}}
train score: 0.93845
test score: 0.9614
time: 1330.1183042526245

たった5epochですが，最適化を行った結果は最適化なしの場合よりかなり良い結果を得ることができていることが分かります．

終わりに

最近は最適化ツールとしてOptunaが人気ですが，様々な最適化手法を試すことのできるツールとして非常に興味深いツールだと思います．
一般にハイパーパラメータ最適化を行うにあたってはOptunaでも使われているTPEが優秀ですが，Meta-Heuristics Algorithmの深層学習のハイパーパラメータ最適化への適用はまだあまり研究が行われていない分野であり，研究を行いやすくするツールとしても役に立つのではないかと思います．