読者です 読者をやめる 読者になる 読者になる

Gunosyデータ分析ブログ

Gunosyで働くデータエンジニアが知見を共有するブログです。

Gunosyデータマイニング研究会 119回, 120回を開催しました

こんにちは。グノシーデータ分析部の関です。 最近はMaison book girlのkarmaをよく聞いています。

今回の投稿では4/24に開催したGunosy DM #119と5/10に開催したGunosy DM #120について紹介します。 これまで同様、これからの強化学習の輪読と論文紹介を行っております。

書籍輪読 - これからの強化学習

#119では2.1節を関が紹介し、 #120では2.2, 2.3節を関が 2.4節をatlimited様に紹介いただきました

1章では価値関数が離散的な状況を想定していましたが、 2.1節では価値関数が連続的であることを考慮し、その中で関数を近似する方法を検討しています。 通常の機械学習では、入力がi.i.dであることを仮定していますが、 強化学習では、得られるデータが方策に依存するので、マルコフ性を持ってしまうため、収束が保証されません。 そのなかで収束をさせるための様々な方法について論じられています。 後半で、すべての手法がセミパラメトリックモデルで解釈できるという点が面白かったです。

2.2節は探索と利用のトレードオフをどのように理解するかという点について紹介されていました。 リグレット、サンプル複雑性などの指標により、評価されています。 ただ、紹介されていたのはシンプルなルールによる手法で、 このあたりはなかなか難しい領域なのかと思うと同時に、 ルールを導入することで、様々な実問題が解決できるヒントでもあるなと思いました。

2.3節は逆強化学習について紹介されていました。 価値関数を定義することが難しい場合における学習法で、 エキスパートや、正しい方策から価値関数を学ぶ方法について学びました。 将棋の棋譜からの学習のように最終的な勝ち負けが与えられた上で、その過程を学習するのは、こういうアプローチなのだなと。

qiita.com

2.4節は経験型強化学習についてでした。 どのように試行回数を減らすかという点についての章で、 学習が過度に行われてしまうようなパターンをうまく回避するような方法が紹介されていました。

論文紹介

#119では 米田より A Complete Recipe for Stochastic Gradient MCMCの紹介が行われました。

この論文はNIPS2015で発表されたもので、一言でいうとすべてのMCMCアルゴリズムが、一般的な確率微分方程式で記述できることを示したものです。 これによって、様々な問題に対するMCMCアルゴリズムを設計することが、容易になったといえます。 当日は数学科出身の米田を中心に、濃密なディスカッションが行われました。

#120では吉田より、 Attention and Engagement-Awareness in the Wild: A Large-Scale Study with Adaptive Notificationsの紹介が行われました。

こちらの論文はYahoo!Japanが発表した論文で、スマートフォンアプリケーションにおけるプッシュ通知の効果を最大化するために、 ユーザの状態を考慮してプッシュを発火させるものです。 実際のYahooアプリにおいて大規模な実験が行われ、その効果が立証されていることが印象的でした。

おわりに

次回の#121は5/16に、#122は5/31に開催予定です。 皆さんの参加をお待ちしております。

Pandasによる実践データ分析入門

こんにちは。データ分析部のオギワラです。最近は「NANIMONO (feat.米津玄師)」をよく聞いています。 今回はPythonのデータ分析ライブラリであるPandasについて、実践的なテクニックを「データ処理」「データ集計(Group By)」「時系列処理」の3カテゴリに分けてご紹介していきます。

Pandasに関する基本的な内容については、前エントリーで既に紹介されているので、是非こちらもご一読して頂けると幸いです。

data.gunosy.io

今回、サンプルデータとしてKaggle上に公開されているFourSquareの行動ログデータを利用しています。

FourSquare - NYC and Tokyo Check-ins | Kaggle

このデータは、2012年の4月から2013年2月まで間にFourSquareを利用したユーザの、以下のようなチェックイン情報が記載されています。

- ユーザID (user_id)
- 入店したお店のカテゴリID (venue_category_id)
- 入店したお店のカテゴリ (venue_category)
- お店の緯度・経度情報 (lat, lng)
- 入店時刻 (timestamp)

このデータを利用して実際の分析内容・シチュエーションにも触れながらPandasの実践テクニックをご紹介していきます。 また、今回使用したJupyter NotebookをGitHubのリポジトリに配置したので、実際の動作等を確認したい方はぜひ利用して下さい。

pandas_analysis_topic/pandas_data_analysis_topic.ipynb at master · ogiogi93/pandas_analysis_topic · GitHub

データ処理

データの取り出し(query)

今回のデータには197ものお店カテゴリが存在し全カテゴリを一度に分析することは困難です。そこで対象のお店カテゴリを交通機関に絞って見ます。そのような時にquery 関数を利用することができます。 query関数はSQLのクエリのように条件を記述することでデータの取り出しを行うことができます。

>>> #venue_categoryがAirtportまたはSubwayである行を抽出する
>>> df.query("venue_category in ['Airport', 'Subway']")
'''
  user_id venue_category  timestamp
10    589 Airport 2012-04-04 04:59:06+09:00
32    1876    Subway  2012-04-04 05:59:52+09:00
34    499 Subway  2012-04-04 06:04:04+09:00
'''

条件文に基づくデータ処理の適用(where)

今回のデータでは特に必要ありませんが、データによって欠損値や異常値が含まれていることがあります。これらの値を含む行を除く場合欠損値はdropna、外れ値はquery関数を利用することができます。しかし行を除くのではなく何かしらの値で置換したいと考えた場合、欠損値は fillna関数で利用することができますが、異常値に対しては利用することができません。このような時にwhere関数を利用することができます。where関数は条件文に基づく行に対してのみ処理を適用することができます。(ex. 負値は0に置換する、異常値は平均値に置換するなど)

>>>#例: user_idが10以下の行に対して, +10加算する
>>>#条件を満たさない行(other)に対して処理が適用されます
>>>df['test'] = df['user_id'].where(cond=lambda x: x >10, other=lambda x: x + 10)
'''
user_idが10以下であるため加算処理が適用されている
  user_id test
9168  1   11
'''
'''
user_idが10以上であるため加算処理が適用されていない
  user_id test
3659  2292    2292
'''

各行への関数の適用(apply)

各行に対して関数を適用したい時、apply関数を利用することができます。また複数の列に対して適用したい時は、行和を指定(axis=1) することで適用することができます。
注意: 大規模なデータセットを分析している場合apply関数は非常に時間がかかるため、DataFrameを一旦リスト化・numpyでベクトル化して関数を適用するなどの工夫をした方が良いです。次回以降の高速化に関する内容も執筆していきたいと思っています。

>>>#user_idをハッシュ化する関数を適用する
>>>df['hashed_user_id'] = df['user_id'].apply(_hash) #user_idを引数とする_hash関数を適用、新たな列を生成

>>> 緯度・経度情報から住所を取得する
>>> df['address'] = df.apply(gis_to_address, axis=1) # 関数上で複数列(緯度・経度)を指定し住所を取得、新たな列を生成

データ集計(Group By)

カラム毎に異なる集計を適用する(agg)

地点・時間別でのチェックイン数の平均・中央値などを知りたい時に、agg関数を利用することでまとめて集計することができます。 今回は例として15分ごとに各店舗(カテゴリ)のチェックイン人数をカウントし、その後店舗(カテゴリ)ごとにチェックイン人数の平均・中央値を求めてみます。

>>># Round関数を適用し、timestampを15分間隔に丸める
>>>df['15min'] = od.DatetimeIndex(df['timestamp']).round('15min')
>>>#15分ごとの各店舗(カテゴリ)のチェックイン人数をカウント
>>> number_of_people = df.groupby(['venue_category', '15min']).agg({'user_id': 'count'}).reset_index()
'''
  venue_category  15min   user_id
10705 Subway  2012-04-04 07:30:00+09:00   8
10706 Subway  2012-04-04 07:45:00+09:00   12
10707 Subway  2012-04-04 08:00:00+09:00   12
10708 Subway  2012-04-04 08:15:00+09:00   12
10709 Subway  2012-04-04 08:30:00+09:00   12
10710 Subway  2012-04-04 08:45:00+09:00   20
'''
>>># 各店舗(カテゴリ)ごとにチェックイン人数の平均・中央値を算出
>>> number_of_people('venue_category').agg({'user_id': ['mean', 'median']})
'''
venue_category    user_id
mean  median
0 Airport 1.207921    1
1 American Restaurant 1.055556    1
2 Antique Shop    1.000000    1
'''

最大・最小値である行を取り出す(first)

各ユーザの最終アクションに関するデータを抽出したい時、groupby及びfirst関数等を用いることで実現することができます。 具体的な流れは(1) 対象のカラム値でソートする(最大値は降順、最小値は昇順) (2) 最大・最小値を抽出したいカラムでグループ化(3) 各グループの最初の行を抽出する 以下に各ユーザの最終地点(timestampが最後の行)を抽出する例を示します。

>>>#(1) 時間順で降順でソート (2) ユーザIDでグループ化 (3) ユーザごとに最初の行を抽出
>>>df.sort_values('timestamp', ascending=False).groupby(['user_id'], as_index=False).first()
'''
  user_id venue_category  timestamp
0 1   Train Station   2012-04-08 12:24:53+09:00
1 2   Train Station   2012-04-13 11:03:36+09:00
2 3   Smoke Shop  2012-04-13 07:26:35+09:00
3 4   Arcade  2012-04-14 11:24:53+09:00
4 6   Train Station   2012-04-08 21:36:10+09:00
'''

標準化や正規化処理を適用する(transform)

カテゴリごとに標準化・正規化したい時、groupby及びtransform関数を利用することができます。 以下に店舗(カテゴリ)ごとにチェックイン数を標準化する例を示します。

>>>#標準化関数
>>>zscore = lambda x: (x - x.mean()) / x.std()
>>>number_of_people['standarized_num_people'] = number_of_people.groupby('venue_category').transform(zscore)
'''
  venue_category  15min   num_people  standarized_num_people
10703 Subway  2012-04-04 07:00:00+09:00   3   -0.501261
10704 Subway  2012-04-04 07:15:00+09:00   4   -0.208489
10705 Subway  2012-04-04 07:30:00+09:00   8   0.962599
'''

時系列処理

時間の丸め処理(round)

既に上記でも登場していましたが、日付や1時間単位での集計だけではなく、15分や3分など独自の時間間隔で集計したいことがあります。このような時、round関数を利用することができます。 round関数はtimestampを指定した間隔に丸めることができます。

>>>df['date'] = df['timestamp'].dt.date
>>>df['hour'] = df['timestamp'].dt.hour
>>>df['15min'] = pd.DatetimeIndex(df['timestamp']).round('15min') #15分間隔
'''
  user_id venue_category  timestamp   date    hour    15min
0 1541    Cosmetics Shop  2012-04-04 03:17:18+09:00   2012-04-04  3   2012-04-04 03:15:00+09:00
1 868 Ramen / Noodle House    2012-04-04 03:22:04+09:00   2012-04-04  3   2012-04-04 03:15:00+09:00
2 114 Convenience Store   2012-04-04 04:12:07+09:00   2012-04-04  4   2012-04-04 04:15:00+09:00
'''

時系列データの差分計算(diff), 変化率計算(pct_change), 移動平均値の計算(rolling)

時系列データの変動を分析したい時、前時刻との差分はdiff関数, 変化率はpct_change, 移動平均値はrolling関数を用いることで簡単に算出することができます。また、店舗(カテゴリ)別やuser_id別で算出したい場合もgroupbyと組み合わせることで実現することができます。

>>>count_people_div_15min['diff'] = count_people_div_15min['num_of_people'].diff() #差分
>>>count_people_div_15min['pct_change'] = count_people_div_15min['num_of_people'].pct_change() #変化率
>>>count_people_div_15min['pct_change'] = count_people_div_15min['num_of_people'].rolling(window=3, center=False).mean() #移動平均
'''
15min num_of_people   change  pct_change  rolling_mean
0 2012-04-04 03:15:00+09:00   2   NaN NaN NaN
1 2012-04-04 04:15:00+09:00   5   3.0 1.5 NaN
2 2012-04-04 04:30:00+09:00   1   -4.0    -0.8    2.666667
3 2012-04-04 04:45:00+09:00   2   1.0 1.0 2.666667
4 2012-04-04 05:00:00+09:00   2   0.0 0.0 1.666667
'''
>>>#店舗カテゴリ別で算出したい場合はgroupbyと組み合わせる
>>>count_people_div_15min_venue_category['change'] = count_people_div_15min_venue_category.groupby('venue_category')['user_id'].diff() #店舗カテゴリ別差分
>>>count_people_div_15min_venue_category['pct_change'] = count_people_div_15min_venue_category.groupby('venue_category')['user_id'].pct_change() #店舗カテゴリ別変化率
>>>_rolling_mean = count_people_div_15min_venue_category.groupby('venue_category')['user_id'].rolling(window=3, center=False).mean() #店舗カテゴリ別移動平均
>>>count_people_div_15min_venue_category['rolling_mean'] = _rolling_mean.reset_index(level=0, drop=True)
'''
venue_category    15min   user_id change  pct_change  rolling_mean
10701 Subway  2012-04-04 06:00:00+09:00   2   NaN NaN NaN
10702 Subway  2012-04-04 06:45:00+09:00   2   0.0 0.000000    NaN
10703 Subway  2012-04-04 07:00:00+09:00   3   1.0 0.500000    2.333333
10704 Subway  2012-04-04 07:15:00+09:00   4   1.0 0.333333    3.000000
10705 Subway  2012-04-04 07:30:00+09:00   8   4.0 1.000000    5.000000
'''

Pandasと合わせて利用したい、おすすめライブラリ

Pandas及びJupyter Notebookでデータ分析していく際に是非利用した方が良いライブラリを2つご紹介します。

for文やapply関数などの進捗をバーで表示する(tqdm)

GitHub - tqdm/tqdm: A fast, extensible progress bar for Python and CLI

Pythonのfor文などの進捗を色付きバーでリアルタイムに表示してくれます。平常心を保ちながら大規模データを分析していくためにはこのライブラリは必須ですね。 Pandas向けの進捗バーも準備されており、tqdmを有効化した後にapply → progress_apply, map → progress_map に変更し実行することで表示されます。

>>> from tqdm import tqdm, tqdm_notebook
>>> tqdm_notebook().pandas() # pandas向けtqdmを有効化
>>> df.progress_apply(f)

f:id:ogiogi93:20170508113353p:plain

DataFrameをMarkdown形式で出力する(pytablewriter)

GitHub - thombashi/pytablewriter: A python library to write a table in various formats: CSV / HTML / JavaScript / JSON / LTSV / Markdown / MediaWiki / Excel / Pandas / Python / reStructuredText / SQLite / TOML / TSV.

分析結果をGithubのissueなどに転記する時に利用しています。DataFrameを自動的にMarkdown形式で出力してくれます。

>>>import pytablewriter
>>>writer = pytablewriter.MarkdownTableWriter()
>>>writer.from_dataframe(count_people_by_15min_category.head(5)) # Markdownに変換したいDataFrameを入力する
>>>writer.write_table()
'''
venue_category|         15min          |user_id|change|pct_change|rolling_mean
--------------|------------------------|------:|-----:|---------:|-----------:
Airport       |2012-04-04T05:00:00+0900|      1|   NaN|       NaN|         NaN
Airport       |2012-04-04T07:00:00+0900|      1|     0|         0|         NaN
Airport       |2012-04-04T08:15:00+0900|      1|     0|         0|        1.00
Airport       |2012-04-04T10:30:00+0900|      1|     0|         0|        1.00
Airport       |2012-04-04T11:45:00+0900|      2|     1|         1|        1.33
'''

さいごに

Pandasには今回紹介したメソッド以外にも、数多くのメソッドが実装されています。ぜひPandasの公式ドキュメントも読んでみて下さい! また、最近Pandasのリポジトリ上に「Pandas Cheat Sheet」という基礎内容が一つのpdfファイルにまとまったシートが公開されているので、ぜひダウンロードしておくことをお勧めします!

github.com

参考資料

【これからの強化学習】 Gunosy データマイニング研究会 #118 を実施しました

gunosy-dm.connpass.com

こんにちは。グノシー開発部のアルシャマンです。最近は、KID FRESINOのSalve feat. JJJをよく聴いています。

今日は4/12(水)に開催したGunosy DM #118について紹介します。前回に引き続きこれからの強化学習の1.3~1.5節の輪読と、論文紹介を行いました。

Gunosy DMこれからの強化学習については、以下のブログ記事で紹介しています。 data.gunosy.io

書籍輪読(これからの強化学習)

データ分析部の大曽根と吉田からそれぞれ1.3~1.4節と1.5節についての発表がありました。

1.3節では、MDP(マルコフ決定過程)における価値関数の表現と、それを推定するアルゴリズムについて学びました。具体的には、ある方策πのもとでの行動価値関数について成立する再帰式であるベルマン方程式とSarsaという学習法、最適行動価値関数について成立する再帰式であるベルマン最適方程式とQ-learningという学習法についてです。

1.4節では、方策を行動価値関数とは別のパラメータで表現された確率モデルと考え、そのパラメータについて最適化することで強化学習問題を解く方策勾配に基づく方法について学びました。

1.5節ではPOMDP(部分観測マルコフ決定過程)における強化学習について学びました。「エージェントは状態を部分的にしか観測できない」とした場合の強化学習についてです。どの状態にいるのかを表す確率分布である信念状態を導入して、価値関数を表現します。

論文紹介

次にデータ分析部の関から、Optimizing the Recency-Relevancy Trade-off in Online News Recommendationsという論文の紹介がありました。

ニュースサイトのトップページに掲載するニュースを選択する時に考慮される「Recency(新しさ)」と「Relevancy(ニュースの重要さ)」というトレードオフ関係にある2つの指標についての分析と「Highest Future-Impact(ニュース掲載後にどれだけPVを稼ぐか)」に基づいたニュース選択方法の提案がされています。

最後に

今後もGunosy DMは隔週開催予定です。次回の#119は、4/25(火)を予定しています。

次回の内容は、今回論文発表を担当した関からこれからの強化学習の2.1~2.3節の輪読と、同じくデータ分析部の【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装までというブログの著者のmathetakeからA Complete Recipe for Stochastic Gradient MCMCという論文の紹介となっています。

興味ある方は、ぜひご参加下さい。

gunosy-dm.connpass.com

【これからの強化学習】 Gunosy データマイニング研究会を実施しました

gunosy-dm.connpass.com

こんにちは、データ分析部の阿部です。 今回は、先日開催したデータマイニング研究会という勉強会についてご紹介します。

データマイニング研究会とは

本勉強会では書籍の輪読と論文紹介を行い、データマイニングに関する基礎知識の向上及び、先端事例の共有・議論を行うことを目的としています。 2週間に1回のペースで開催されており、社外にも公開し広く知見を共有することを目指しています。 Gunosy創業時から取り組んでいるためこの手の勉強会としては歴史は長く(?)、今回で117回目になりました。

f:id:y-abe-hep:20170407102639j:plain

f:id:y-abe-hep:20170407102646j:plain

これからの強化学習

今回からは「これからの強化学習」を進めていて、1.1と1.2を終わらせました。 内容は強化学習の基礎的なところで、強化学習の構成要素が中心となっています。

これからの強化学習

これからの強化学習

  • 作者: 牧野貴樹,澁谷長史,白川真一,浅田稔,麻生英樹,荒井幸代,飯間等,伊藤真,大倉和博,黒江康明,杉本徳和,坪井祐太,銅谷賢治,前田新一,松井藤五郎,南泰浩,宮崎和光,目黒豊美,森村哲郎,森本淳,保田俊行,吉本潤一郎
  • 出版社/メーカー: 森北出版
  • 発売日: 2016/10/27
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログ (2件) を見る

本書はこれから研究を始める大学院生や研究者に向けて、最新の研究を理解する手がかりになることを目指しているとあります。 実務でも強化学習は応用されており、例えば広告のクリエイティブ選択に強化学習の一種とも言えるバンディットアルゴリズムが使われていたりします。 この他にも応用できそうな箇所もあるので、将来的に本書の内容を活かしてサービスへの応用もできればと考えています。

ちなみに、前回までは「トピックモデルによる統計的潜在意味解析」を進めていました。 こちらも良書なので、紹介しておきます。

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

論文紹介

論文紹介はこちらです。

Focus on the Long-Term: It’s better for Users and Business(KDD2015)

Focus on the Long-Term: It's better for Users and Business

Web系企業ではA/Bテストによってプロダクトの品質を分析するケースが多いですが、どうしても短期的な指標の評価に重点を置いてしまい長期的な価値が評価しきれていないこともあります。 この論文ではどのようにして長期的価値を評価するのかということを論じており、興味深かったです。

最後に

今回はGunosyデータマイニング研究会について紹介しました。 興味があったらぜひご参加ください!

【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装まで

こんにちは。初めまして。
データ分析部新入りのmathetake(@mathetake)と申します。

先日個人ブログでこんなエントリを書いた人です:

mathetake.hatenablog.com

そんなこんなでTwitter就活芸人(?)として活動(?)してましたが、これからは真面目に頑張っていこうと思います。



今日はみんな大好きベイズモデリングおいて、事後分布推定に欠かせないアルゴリズム(群)の一つである*1

マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo)

通称MCMCに関するエントリです。より具体的に、

MCMCの意義(§1.)から始め、マルコフ連鎖の数学的な基礎(§2.,3.,4.)、MCMCの代表的なアルゴリズムであるMetropolis-Hastings法(§5.)、その例の1つである*2Langevin Dynamics(§6.)、そして(僕の中で)絶賛大流行中のライブラリEdwardを使ってより発展的(?)なアルゴリズムであるStochastic Gradient Langevin Dynamicsの説明&実装(§7.,8.)していきたいと思います。


今までのデータ分析ブログのエントリと少しテイストが違うかもしれませんが、お楽しみいただけたら幸いです。

細心の注意を払ってはいますが、数学的正しさを保証する記事ではありませんので、詳細が気になる方は§Referenceにある資料を御覧ください。

はじめに

本題に入る前にまず、MCMCの意義について軽く触れておきます。

ベイズ統計におけるモデルパラメータの推定には事後分布からのサンプリングがかかせません。即ち、モデルのパラメータを{\theta \in \Theta},データの集合を{\mathcal{D}}として、ベイズの定理から得られる次のような確率分布

{ p \left( \theta \ |  \ \mathcal{D} \right) = \dfrac{p \left( \mathcal{D} \ | \ \theta  \right)p(\theta) }{ \displaystyle \int_\Theta p \left( \mathcal{D} \ | \ \theta  \right) p(\theta) d\theta } }

からサンプリングする必要がありますが、一般に分母( 正規化定数 )は解析的に求められず、サンプリングは困難です。

そこを、マルコフ連鎖の収束定理等の数学的基礎理論を使って上手くサンプリングする手法がMCMCで、広く応用されています。*3

ベイズ統計の数学的背景については、僕の個人ブログのエントリを御覧ください

mathetake.hatenablog.com

離散値マルコフ連鎖

{ \mathcal{X}:=\{X^0,X^1,X^2,\dots,X^t,\dots \} } を離散な集合{ \{1,2,3,4,..,k\}(k \in \mathbb{N}) } に値を取る、離散な確率変数の集合とします。確率的な値の時間発展だと考え、確率過程と呼ぶことにします。

(Definition 1.) {\mathcal{X} }マルコフ連鎖であるとは、全ての{ \{x^{s}\}_{s=0}^{t+1} \subset \{1,2,\dots,k \}}に対して

{ \mathbb{P} \left(X^{t+1} = x^{t+1} \ |  \ X^t=x^t ,\dots X^0=x^0 \right) = \mathbb{P}\left(X^{t+1} =x  \ | \ X^t=x^t  \right)\ \ \   \cdots (1) }

が成立する事。

感覚的には、次の時刻における値の分布は現在の値のみで決まり、それ以前の値には影響されない確率過程の事です。

数学的便宜上、次の定義も用意しておきます。

(Definition 2.) マルコフ連鎖{\mathcal{X} }斉時的であるとは、

{p_{i,j}^t:= \mathbb{P}\left(X^{t+1} =j \ | \ X^{t} = i \right) }

時間に依存しない事。またこの時 行列{ p_{i,j} := p_{i,j}^t }遷移行列 と呼ぶ。

つまり状態遷移の確率が時間に依らず一定なマルコフ連鎖の事で、状態遷移を表す行列を遷移行列 としています。


ずっと数学の話になってしまいましたが、最後に一つだけMCMCに関わる重要な命題とその系を述べておきます。

(Proposition 3.) {\mathcal{X} }を斉時的なマルコフ連鎖とする。この時 {X^t} の分布は 遷移行列{ p_{i,j} }{X^0}の分布(初期分布)により完全に決定される。

(Proof) {X^0,X^1,\dots} の分布を{\pi^0, \pi^1,\dots,}とし、{\pi^t_j}で確率変数{X^t}が値{j} を取る確率とする。マルコフ連鎖の性質(1)より

{ \begin{eqnarray} \pi^1_j &=& \sum_{i=1}^k \mathbb{P} \left( X^1=j \ | \ X^0=i \right) \mathbb{P} \left(  X^0 =i \right)  \\
&=& \sum_{i=1}^k \pi^0_i p_{i,j} \ \ \  \cdots (2) \end{eqnarray}}

また、

{ \begin{eqnarray} \pi^t_j &=& \sum_{i=1}^k \mathbb{P} \left( X^t=j \ | \ X^{t-1}=i \right) \mathbb{P} \left(  X^{t-1} =i \right)  \\
&=& \sum_{i=1}^k \pi^{t-1}_i p_{i,j}  \ \ \  \cdots (3) \end{eqnarray} }

が成立するので、帰納的に示される □

この命題の系として次が得られます

(Corollary 4.) 初期分布{\pi^{0}} と遷移行列 {T:= (p_{i,j})} を与えることで、式(2),(3)により斉時的マルコフ連鎖が得られる。

マルコフ連鎖収束定理とMCMC

あと少しMCMCの説明まで辿りつきます。もう少々数学にお付き合いください。

§1. で紹介したように、サンプリングしたい確率分布{ \pi =(\pi_i) }が手元(?)にあるとします。

(Definition 5.)遷移行列が{T:= (p_{i,j})} により与えられる、斉時的マルコフ連鎖{ \mathcal{X}}{\pi}不変分布に持つとは

{ \displaystyle \sum_{i=1}^k \pi_i p_{i,j} = \pi_j  \ \ \ ( \ \forall j \in \{ 1,\dots,k\})}

が成立する事。行列の式で書けば

{ \pi T = \pi  }

が成立すること。この時{\pi}{ \mathcal{X}}不変分布と言う。

さて、MCMCの肝となる定理は次のものです

(Theorem 6.(離散値マルコフ連鎖の収束定理))
{\pi } を不変分布に持つ斉時的マルコフ連鎖{ \mathcal{X}}(i) 非周期的 かつ (ii) 既約*4 である時、マルコフ連鎖は不変分布 {\pi } に収束する。即ち


{ \displaystyle \lim_{t \rightarrow \infty} \mathbb{P}\left( X^t = i \right)  =  \pi_i  }

が成立する。

ここまで来てやっと、MCMCの定義を与えることができます:

(Definition 7.) MCMC(Markov chain Monte Carlo)とは、サンプリングしたい確率分布{\pi} を不変分布とするような既約で非周期的なマルコフ連鎖を構築&サンプリング するアルゴリズムの事。

Theorem 6.により、MCMCにより生成されるサンプルの列{ \{x^t\}^\infty_{t=0} } は確率分布 {\pi}からのサンプルに収束し、目的を達成することができます。

既約性と非周期性を満たすようなマルコフ連鎖を構築するのはそんなに難しくはありません、が、サンプリングしたい確率分布{\pi} を不変分布とするようなマルコフ連鎖を構築するのは一般に困難です。

そこでよく用いられるのが*5詳細釣り合い条件と呼ばれる不変分布を持つための十分条件です:

(Definition 8.)遷移行列が{T:= (p_{i,j})} により与えられる、斉時的マルコフ連鎖{ \mathcal{X}}不変分布 {\pi} に対して詳細釣り合い条件を満たすとは


{ \displaystyle  \pi_i \  p_{i,j} = \pi_j \ p_{j,i}  \ \ \ \ (  \forall \ i, j \in \{ 1,\dots,k\})}


が成立する事。またこの時、マルコフ連鎖{ \mathcal{X}}{\pi} を不変分布に持つ。

連続値マルコフ連鎖の場合

今まで簡単のため、離散な値を持つマルコフ連鎖について話をしてきましたが、連続な確率変数の話に一般化する事ができます。

連続値( 便宜上{\mathbb{R}^d}値とする )なマルコフ連鎖{\mathcal{X} = \{ X^0,X^1,\dots,X^t,\dots \}}斉時的である時、離散値の場合の推移行列に対応する推移核{ T: \mathbb{R}^d \times \mathbb{R}^d \rightarrow \mathbb{R} }

{  \mathbb{P} \left( X^{t} \in A \  | \ X^{t-1} =x \right) = \displaystyle \int_A T(x,y) dy \ \ \ \ \ \ (A \subset \mathbb{R}^d ) }

を満たすものとして定めます。感覚的には現時刻で値 {x} を取る時、次の時刻の分布を表す密度関数です。

また、{ \mathcal{X}}が分布{\pi}を不変分布に持つとは

{ \displaystyle \int_{\mathbb{R}^d} \pi(x)T(x,y)dx = \pi(y) }

を満たすことであり、詳細釣り合い条件

{ \pi(x) \  T(x,y) = \pi(y) \ T(y,x)   \ \ \ \cdots (4)}

で与えられます。

(Theorem 9.(連続値マルコフ連鎖の収束定理))
{\pi } を不変分布に持つ斉時的マルコフ連鎖{ \mathcal{X}}(i) 非周期的 かつ (ii) 既約である時、{ \mathcal{X}}total variation distanceの意味で{\pi} に収束する。即ち、


{ \displaystyle \lim_{t \rightarrow \infty} \sup_{A \subset \mathbb{R}^d } \left| \pi^t(A) - \pi(A) \right| \ = \ 0 }


が成立する。

Metropolis-Hastings法

MCMCの代表的なアルゴリズム(群)である、Metropolis-Hastings法(以下M-H法)について説明します。

M-H法では、各 {x}に対して提案分布と呼ばれる確率分布 {q(x,y)dy}を用意し、採択確率と呼ばれる確率

{\alpha(x,y):= \min \left\{ 1, \dfrac{\pi(y)q(y,x) }{\pi(x)q(x,y)} \right\} }

を準備します。そして推移核 {T(x,y)}

{T(x,y) = \alpha(x,y)q(x,y) + A(x)\delta_{x,y} }

として定義&適当な初期分布{\pi^0 }を与えることで斉時的マルコフ連鎖 {\mathcal{X}_{MH} }を考えます。(ここで {\delta_{x,y} } はクロネッカーのデルタ関数で、{A(x)} は正規化定数を与える関数。)

(Theorem 10.) 斉時的マルコフ連鎖 {\mathcal{X}_{MH} }{\pi} を定常分布に持つ

(Proof) 定義から

{ \begin{eqnarray} 
\pi(x)q(x,y)\alpha(x,y) &=& \pi(x)q(x,y) \min \left\{ 1, \dfrac{\pi(y)q(y,x) }{\pi(x)q(x,y)} \right\}\\
&=& \min \left\{ \pi(x)q(x,y), \pi(y)q(y,x) \right\}\\
&=& \pi(y)q(y,x) \min \left\{ \dfrac{\pi(x)q(x,y)}{\pi(y)q(y,x)}, 1 \right\}\\
&=& \pi(y)q(y,x)\alpha(y,x)
\end{eqnarray}}

が従うので、クロネッカーのデルタの定義から(4)式が成立し詳細釣り合い条件が満たされる。□

Theorem 10.だけでは保証されない既約性や非周期性を満たすような推移核の具体的な設計は重要な課題ですが、ここではそのような性質が満たされ Theorem 6.(or 7.) が成立すると仮定しましょう。*6

この時マルコフ連鎖{\mathcal{X}_{MH}} からサンプリングする次のようなアルゴリズムを Metropolis-Hastings法と呼びます。

(Metropolis-Hastings法)
(1) 初期分布{ \pi^0} から {x^0} をサンプリングする:

{x^0 \sim \pi^0 }

(2) {t =0,1,2,3, \dots} に対して、以下を実行する
 (i) 標準一様分布から乱数 {u} を生成する:

{u \sim U(0,1) }
 (ii) {y} を提案分布 {q(x^t, y)} からサンプリングする:
{y \sim q(x^t, y) }
 (ii) 次の式により"次の点" {x^{t+1} }を決める:
{x^{t+1} \ := \ \displaystyle
  \begin{cases}
    y  \ \ \ \ \ if \ \  u < \alpha(x^t, y) \\
    x^t \ \ \ \ otherwise
  \end{cases}
}

このアルゴリズムが、実際に上で与えた{\mathcal{X}_{MH} } からサンプリングしている事は明らかでしょう。

§0. で述べたように採択確率の計算には目標の分布{\pi(x))}正規化定数が必要ない事に注目して下さい。

Langevin Dynamics

M-H法を実際に実行するためには、提案分布 {q(x,y)dy} を定義する必要があります。

ここではその例として、Langevin Dynamics法(Metropolis-adjusted Langevin Algorithm)(以下LD法) を紹介します。

LD法では提案分布を次のように定義します:

{q(x,y) := \mathcal{N} \left( x + \dfrac{\epsilon}{2} \dfrac{\partial \log \pi }{\partial x}(x),  \ \ \epsilon I \right)  }

ここで{\mathcal{N}\left( \mu, \Sigma \right)} は平均 {\mu}, 分散 {\Sigma}に従う正規分布で、{\epsilon}はstep size(またはlearning rate)と呼ばれるハイパーパラメータです。

ベイズモデルの事後分布に適応する場合において、決定論的に眺めると、LD法は事前分布で正規化した、ノイズ入り勾配降下法のようなものであると解釈することができます。

実際、事後分布推定において、{\pi} は観測データを {\mathcal{D = \{d_i \}}_{i=1}^N  }として

{ \pi(\theta) = p \left( \theta \ |  \ \mathcal{D} \right) = \dfrac{p(\theta) \prod_{i=1}^N p\left( d_i \ | \ \theta  \right) }{ \displaystyle \int_\Theta p \left( \mathcal{D} \ | \ \theta  \right) p(\theta) d\theta } }

で与えられ、その対数微分は

{\displaystyle \dfrac{\partial}{\partial \theta} \log \pi \left( \theta \right) =\dfrac{\partial}{\partial \theta} \log p \left( \theta \right)  + \dfrac{\partial}{\partial \theta} \sum_{i=1}^N \log p\left(d_i \ | \ \theta \right)  \ \ \ \cdots (5)}

のように計算できるので、上のような解釈ができます。

実はLG法は、Stanで有名になったHamiltonian Monte Carlo法の特別なケースと等価になっているので、気になる方は [1]や[2] を御覧ください。

Stochastic Gradient Langevin Dynamics(SGLD)

ここまでMCMCの数学的基礎からM-H法、そしてその具体例としてLG法を紹介しました。

LG法の問題点として、(5) 式の計算量がサンプル数が増えるほど膨大になっていく点があります。

近年はビッグデータと呼ばれるバズワードもあるように、サンプル数&パラメータ数が巨大なセッティングでモデリングする事が多いのでこのままLG法を適用する事はできません。*7


その問題点を克服するサンプリング手法として、ここで紹介するのが Stochastic Gradient Langevin Dynamics法 [8](以下SGLD法)です。

SGLD法では次のように初期分布からサンプルしたパラメータを更新していきます:

まず、次の2つの条件

{\displaystyle \sum_{t=0}^\infty \epsilon_t^2  < \infty \ , \ \ \ \ \ \ \sum_{t=0}^\infty \epsilon_t  = \infty}

を満たす数列{ \{\epsilon_t \}_{t=0}^\infty }を用意し、パラメータサンプル{ \{ \theta_t \}_{t=0}^\infty }を次の式によって取得して行きます

{\begin{eqnarray}
\displaystyle & \theta_{t+1} &\leftarrow  \ \ \ \theta_t  +  \dfrac{\epsilon_t}{2} \dfrac{\partial L}{\partial \theta} (\theta_t) + \eta_t  \ , \ \ \ \eta_t  \sim \mathcal{N}\left( 0 \ ,\epsilon_t \right) \\ \ \\
&\dfrac{\partial L}{\partial \theta} (\theta_t)& =  \dfrac{\partial \log p}{\partial \theta_t}  \left( \theta_t \right) +  \dfrac{N}{\left| \mathcal{S_t} \right|}\sum_{d \in \mathcal{S_t} } \dfrac{\partial \log p \left(d \ | \ \theta \right) }{\partial \theta}(\theta_t)
\end{eqnarray} }

ここで、{\mathcal{S_t}} はデータ{ \{d_i\}_{i=1}^N } からランダムに抽出された {N} より十分小さいミニバッチとします。

注意として、このアルゴリズムに対応するマルコフ連鎖は斉時的ではないので、上述の収束定理は適用できません。

ですが、例えば [9]で収束性に関する解析がされています。

EdwardでのSGLDの実装

最後に確率モデリング用ライブラリEdwardを用いて、SGLDをベイズ的線形回帰に適用してみようと思います。

Edwardの詳しい使い方は公式チュートリアルまたは次の論文

[1701.03757] Deep Probabilistic Programming

[1610.09787] Edward: A library for probabilistic modeling, inference, and criticism

をご覧ください。また質問等ありましたら@mathetakeまで気軽にリプライorDMください。

まず各種ライブラリをimportします。

import numpy as np
import tensorflow as tf
import edward as ed
from edward.models import Normal, Empirical
import time

次にデータセットを用意します。

N = 20000  # サンプル数
D = 50  # 特徴量の次元
N_ITER = 10000  # MCMCのiteration
MINI_BATCH_SIZE = 2500  #ミニバッチのサイズ 

# toy dataset. 切片=0はなし.
def build_toy_dataset(N, D, noise_std=0.1):
    w = np.random.randn(D).astype(np.float32)
    X = np.random.randn(N, D).astype(np.float32)
    Y = np.dot(X, w) + np.random.normal(0, noise_std, size=N)
    return w, X, Y

# データ生成。観測値のノイズの分散は既知とする。
w_true, X_data, Y_data = build_toy_dataset(N, D)

# ミニバッチを返す関数 
def next_batch(mini_batch_size=128):
    indexes = np.random.randint(N, size=mini_batch_size)
    return X_data[indexes], Y_data[indexes]

モデルを構築し, 推論のためのインスタンスを作ります。

# 観測データを挿入するためのデータを収めるplaceholder
x = tf.placeholder(tf.float32, [MINI_BATCH_SIZE, D])
y_ph = tf.placeholder(tf.float32, [MINI_BATCH_SIZE])

w = Normal(mu=tf.zeros(D), sigma=tf.ones(D))
b = Normal(mu=tf.zeros(1), sigma=tf.ones(1))
y = Normal(mu=ed.dot(x, w) + b, sigma=tf.ones(MINI_BATCH_SIZE)*0.1)

# 経験分布をposteriorの近似に使う
qw = Empirical(params=tf.Variable(tf.random_normal([N_ITER, D])))
qb = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 1])))

# SGLD法用インスタンス
SGLD = ed.SGLD(latent_vars={w: qw, b: qb}, data={y: y_ph})


最後に推論を実行します。

# 推論GO 
# data辞書にはobservedな確率変数の観測データを送る。
# xの値は確率変数ではないので、updateの際feed_dictで送る。
SGLD = ed.SGLD(latent_vars={w: qw, b: qb}, data={y: y_ph})
SGLD.initialize(scale={y: float(N) / MINI_BATCH_SIZE}, step_size=0.00001, n_iter=N_ITER)

start = time.time()
init = tf.global_variables_initializer()
init.run()
for _ in tqdm(range(N_ITER)):
    X_batch, Y_batch = next_batch(MINI_BATCH_SIZE)
    _ = SGLD.update(feed_dict={x: X_batch, y_ph: Y_batch})
elapsed_time = time.time() - start
print("elapsed_time:{}".format(elapsed_time))


実行結果ですが、ミニバッチで勾配計算をしない通常のLD法と比較してみました:

アルゴリズム 実行時間(s) MSE of W Estimated b
SGLD法 26.4 {4.9 \times 10^{-5}} 0.0020
LD法 55.4 {5.3 \times 10^{-7}} 0.0001

といった感じです。精度はLD法には劣りますが、実行時間は約2倍短いと言った感じです。

より大きなデータセット&複雑なモデルの場合、この差はより顕著になるでしょう。

Appendix

(Definition) 斉時的マルコフ連鎖{ \mathcal{X}}非周期的であるとは、任意の{i \in \{1,\dots,k\} } に対して

{gcd \left\{ n \ | \ \mathbb{P} \left( X^n=i \  |  \ X^0 =i \right)  > 0 \right\} = 1 }

が成立する事。

(Definition) 斉時的マルコフ連鎖{ \mathcal{X}}既約であるとは, 任意の{i, j \in \{1,\dots,k \}} に対してある{t_{i,j} \in \mathbb{N} } が存在して

{ \mathbb{P}\left(X^{t_{i,j}}=j \ | \ X^0 =i  \right) >0   }

が成立する事。

*1:他には変分ベイズ法などがあります

*2:Hamiltonian Monte Carlo法の一種でもあります。

*3: 事後分布推定に限らず、”正規化定数が分からない分布からサンプリングする手法” として広く使われています。

*4: (i) 非周期的 かつ (ii) 既約 の定義はAppendixに付けておきます。

*5:詳細釣り合い満たさないようなMCMCの研究が最近流行っている、らしいです。例えば 詳細つりあいを満たさないマルコフ連鎖モンテカルロ法とその一般化を御覧ください。

*6:M-Hについて既約性や非周期性が満たされるための条件は論文[5]をご覧ください。

*7:LG法にかぎらず、M-H法は採択確率の計算が基本的にintractableです。ですが、例えばハミルトニアンモンテカルロ法をビッグデータに適用する手法については次の論文があります : Stochastic Gradient Hamiltonian Monte Carlo