2017-11-03

アルゴリズムレベルの並列分散処理の設計

はじめに

私は約3年間並列分散処理を専攻してきました。データマイニング系の研究室であったため、周囲に同様の分野を専攻する人がおらず、基礎を確立させるには辛いものがありました。並列分散処理に対する自分の知見を周りと共有したことがあまりないので、本当に自分の持つ知識は正しいのか不安です。ここで、3年間の知識を簡単に整理することにしますが、話の内容は多岐にわたっているので文章にあまりまとまりがないかもしれません。主にアルゴリズムの並列化・並列分散化についてまとめました。なお、私は並列分散処理基盤に関して Hadoop と Spark しか触れておらず、主に OpenMP のような低層の話はよく理解していません。

並列処理と並列分散処理の基礎

言葉の定義

まず並列処理とは、1台の計算機が持つ複数のコアを用いて処理を並列に実行すること。分散処理とは、複数台のサーバーを用いて処理を実行すること。並列分散処理はその両方を兼ね備えた処理のことを指すとする。Intelなどのハイパースレッディング・テクノロジーを搭載している場合は、一つのコアを複数のコアと認識できるため、並列に処理できる数を増やせる。ただし、その性能は複数コア分にはならないようなので、厳密な扱いには注意が必要。

並列化と分散化

並列化と分散化は似たようなところはあるかもしれません。直列処理から並列処理への移行は、一つのコアで処理させていたものを複数のコアで処理させる、といったように概念の変化があります。一方で、直列処理もしくは並列処理を分散化させるには、一つの計算機で処理したものを複数の計算機で処理する、といったように概念の変化があります。まず分散処理は一つのタスクを複数のタスクに分割して問題を解くため、分散処理を取り組む前に並列処理を取り組む必要があるでしょう。ただ、並列処理はコア単位、分散処理は計算機単位の話になるため、二つの特性は大きく異なります。

並列処理の背景

並列処理は複数のコアがあって、初めて効率化ができます。現在は、一般のユーザーが8コアのノートPCを使用する時代ですが、このようなマルチコアはどういうときに有効に働くのだろう。例えば、二つのアプリケーションを起動していたならば、コア1でアプリケーション X を実行し、コア2でアプリケーション Y を実行できます。このとき、並列の粒度はアプリケーション単位と言えるでしょう。ソフトウェアエンジニアは並列化という概念を全く知らずにプログラムを実装したとしても、ユーザーはアプリケーションを複数立ち上げることはしばしばあるので、その恩恵を大いに受けるでしょう。

直列実装の効率性の低さ

初学者が初めて作成するようなライブラリを使わないC言語やJavaのプログラムは、シングルスレッドで動作するプログラムです。例えば、8コア搭載のPCでシングルスレッドのプログラムを動作させる場合には、最大12.5%の性能しか得られないでしょう。まだプログラムの計算量が少ないから気にすることはありません。ただし、問題はとんでもなく重い処理を一つのプログラムかつシングルスレッドで実行している場合です。プログラムの実行から出力を得られるまで数分から数時間かかるようなプログラムをただ待つのは疲れるでしょう。何のライブラリも使用せずに、並列化を考慮しない純粋な実装をしたときには、間違いなくシングルスレッドで動作するプログラムのはずです。並列化を検討しないと、手元にある高い性能を持ったマシンが、マルチコアの恩恵を一切受けることなくのんびりと処理をし続ける羽目になります。現代のコンピュータアーキテクチャーに対して、直列実装は適さないでしょう。

並列に詳しくないとマルチコアの恩恵を受けられないのか

実はそんなこともないようです。プログラムの実装に使用するライブラリが並列処理をサポートしているかもしれません。もしパフォーマンスを改善したいならば、どんなライブラリが並列処理をしているのかは調べるとよいでしょう。使用するライブラリによっては極端な実行時間の差が出るかもしれない。特に低いレイヤーで実装をする人は気をつけたほうがよいのかもしれない。

並列分散化の実行容易性

現在、並列分散処理はライブラリを用いれば簡単に実行できます。これは、Apache Spark のGithubの例に掲載されているPythonのコードを少しだけ変化させたものです。

sc.parallelize(range(1000),5).count()

このコードの意味は、1から1000までの離散値集合を5つのデータに分割し、個々のプロセッサでそのデータを数え上げた後に集約をします。並列分散処理を実行するための基盤があれば、細かい並列分散処理はバックエンドで実行してくれるので、ユーザー側の実装は非常に簡単になります。

じゃあ並列アルゴリズムで問題を解こう

全てのアルゴリズムが並列化できるわけではありません。並列化するにはまず問題を解く計算をどのように分割するのかを考えなければいけません。

並列化可能な問題とは

並列計算の粒度を決めるということです。並列処理の入門がワードカウントです。並列分散処理のフレームワークであるMapReduceなどでも問題の解き方としてよく挙げられる例です。並列計算の粒度を単語とし、それを幾つかのプロセッサーに振り分ける。単語の頻度テーブルを各プロセッサーで作成し、最後にそのテーブルを合成する。このワードカウントの問題を一般化すると、「データをいくつかのチャンクに分割し、そのチャンクごとに問題を解く順番が関係性がないという条件のもとで、非同期的に各プロセッサで処理実行し、その計算結果を集約できる」ということになります。このように問題を定式化できないときに、並列化はできません。

並列化不可能な問題とは

前の計算結果を利用して次の計算をしなければならないような、順次的に結果を更新しなければいけないタイプの問題は並列化不可能です。この場合、f(t)を計算するためにtを必要とし、それが繰り返しで表現されています。

int t = 0;
for(i = 0; i < 10000; i++){
 t += f(t);
}

この例はこれ以上考えようがないように思います。ただし、問題を並列化で解く時に考えるのは、どの部分を並列処理の単位とするかです。異なる部分を対象とすることで並列化できるかもしれません。なので、アルゴリズムの根本を考え直すと、並列化へ持ち込むことも可能な場合があります。

並列分散化を考える前に

素直に問題を解いて処理が重たすぎると感じた場合、計算量が中程度でありその計算を用いて大量のデータを処理しなくてはならない場合、もしくは計算量が多いが有限時間内に計算を解くことができると理解している場合に、初めて並列分散化を検討するべきでしょう。特に問題となっていないのに、むやみに並列分散化から検討するのはよくないと感じます。もしかしたら、純粋に並列化をするだけで問題を解ける可能性があるからです。並列分散化すると早くなるというモチベーションで挑むと厄介なことになります。

並列分散処理で考えるべきこと

アルゴリズムが正しい設計なのかをまず確認するべきです。設計が正しくないのならば、複数の計算機を用いて並列分散処理を実行させることは資源の無駄遣いとなってしまいます。

並列化と分散化の困難性

並列化の問題に限定して考えると、時間計算量や空間計算量などを事前に見積もり、どこがボトルネックになりそうかを判断する知識、直列アルゴリズムを並列アルゴリズムに落とし込むための手法の考案、それを解くためにどういった性能を持つハードウェアを用意するべきなのか、対費用効果で良い結果を得ることができるのか。分散化の観点ではもっと難しくなります。ネットワーク周辺やセキュリティの問題、複数の計算機での効率の良いデータの持ち方、デバッグの考え方といったように、単体の計算機での処理が複数台になるときに非常に大きなパラダイムの変化があります。

考えることが多すぎるが...

CやC++などの低層レベルでプログラムを実装しようとなると、それ以上の困難性が付きまといます。ただ、ここまで検討するのが予測できなくても、ビッグデータというバズワードで流行した並列分散処理の分野では、おおよその知見と事例が積み重ねられてきました。どういう場合に並列分散化を検討するべきなのかは、実例をみることがよいでしょう。技術的な観点からはSpark、AWS Elastic MapReduceやCloudera Manegerといった並列分散処理機構を駆使し、効率よく問題を解くことで、そのユーザーが向き合うべき問題の本質により時間をかけることが可能になってきています。

並列分散化だけがパフォーマンスを改善するわけではない

並列分散化だけがパフォーマンスを改善させるわけではありません。

設計を考える順番

直列化、並列化、並列分散化の順番で設計を考えるべきだと思っています。この順番で規模が大きくなり、それが大きくなるほど、あらゆる側面を考慮して設計をしなければならないからです。

直列実装：並列化困難（並列処理の粒度が細かすぎる）もしくは、小規模の時間計算量と空間計算量を要する場合
並列実装：並列化が可能（並列処理の粒度を大きくできる）かつ、直列実装で処理困難である場合
並列分散実装：並列実装でも処理しきれないレベルの時間計算量と空間計算量を要する場合

ちなみに時間計算量が指数的に増加する場合、並列分散では線形にしか対応できないので、データ量が十分に大きい場合対応不可能です。なので、計算量が指数的に増加する場合、データ量が小さいという条件が必要になるでしょう。

直列化だけで改善することを試みる

実行速度を改善したいならば、より高速に動作するCやC++といったプログラミング言語の選択を検討するのもよいでしょう。実装言語によっても実行時間に100倍スケールの差が出ることもあるからです。感覚的に低層のプログラミング言語は並列化の実装が困難になります。開発効率という観点を評価しなければならず、実装言語を変更できないのであれば、並列化を検討してみましょう。

並列化しても高速化されるとは限らない

アルゴリズムを並列化すると、処理時間がコア数に反比例するというわけではない。アムダールの法則というものがある。99%の並列化であっても、プロセッサ数が99の場合、残りの 50%の性能向上を無駄にします。多くのプロセッサを使用する場合は、基本的には全てのアルゴリズムの並列化を問われます。並列化しても使用したプロセッサ数だけの高速化は見込めないことはまず当然の話になってきます。

まとめ

アルゴリズムの並列分散処理をする場合には、山ほど考えるべきことがあります。並列化を検討してもなお処理が終わらないといった問題に直面した時に、はじめて並列分散処理で実装することを考えるべきでしょう。まだ語りきれない部分も多くありますが、一旦締めとします。

2017-09-30

どの犯罪を犯したのかを分類する（Kaggle体験記）

はじめに

データ分析初心者であるので、kaggleを通して学んだデータ解析の手法をまとめるためにかきます。初学者にもわかるようになるべく丁寧に説明します。

犯罪の分類

Kaggleとは、データ分析のコンペティションサイトです。多くのデータ分析者がデータに基づいてモデルを生成し、未知を予測をします。kaggleでよく話題にされる初学者向けの分析は、タイタニック号の乗船者情報による生存可否の予測や人が書いた数字画像データから数字を分類することなどです。今回は、San Fransico Crime Classification（サンフランシスコにおける犯罪者の罪状の分類）を取り組みます。

データセットを見る

学習データの大きさは、22.09MBで約87万のレコード数になります。

$ wc -l train.csv
878050 train.csv

"Dates"（犯罪の発生日）、"Category"（罪状）、"Descript"（罪状の詳細）、"DayOfWeek"（曜日）、"PdDistrict"（警察署の名前）、"Resolution"（事件の解決方法）、"Address"（犯罪が発生した通り）、"X"（経度）、"Y"（緯度）の9種類の変数で構成されています。その中の"Category"は、今回予測すべき値です。この課題のスコアの基準は、各"Category"に割り当てられて正解の確率値と予測した確率値の対数誤差で計算されています。そのため、例えば、ある犯罪者の"Category"は放火が0.12の確率で、続いて麻薬取引が0.05の確率で...といった形で出力しなければいけません。

犯罪情報を分類をするには？

訓練データとその正解ラベルが事前に与えられているため、この問題を教師あり学習で解きます。今回は、この問題を解く手法としてランダムフォレストを採用します。数学的に難しい実装部分はライブラリに依存するため、考える必要はありません。分類器に入れるデータは数値データでなければいけないため、特徴量をどのようにして処理するかがポイントです。また、今回はkaggleに置かれたカーネル *1を参考にしたため、よく似た形になってしまいましたが、随所調べてまとめました。この後に実行する手順を簡単にまとめると、

データを前処理
分類器のパラメータを指定（各分類器の設定）
データを学習（fitメソッド）
予測（predictメソッド）

となる。

データセットの処理方法を考える

分類器に入れるデータを数値データに変換していくことを念頭に入れながら、データの前処理を考えていきます。

Dates	犯罪の発生日	"yyyy/MM/dd hh:mm:ss"のフォーマットを年・月・日・時のカテゴリカルデータに変換
Category	罪状	目的変数であるため除外
Descript	罪状の詳細	内容が細かすぎるため、使用不可
DayOfWeek	曜日	日月...土、とカテゴリカルデータに変換
PdDistrict	警察署の名前	種類がかなり限られているためカテゴリカルデータに変換
Resolution	事件の解決方法	Null値が多すぎるため、使用不可
Address	犯罪が発生した場所	種類がかなり限られているためカテゴリカルデータに変換
X	経度	数値データとしてそのまま使用
Y	緯度	数値データとしてそのまま使用

まず、"Descript"は、記述内容が多岐にわたるため、自然言語処理を含めないと処理ができないと考え、このデータを使用しないことにしました。"Resolution"は、Noneの値が多すぎるため、このデータを使用することは難しいと判断します。"Dates"と"DayOfWeek"は、カテゴリカルデータとして扱うため、前処理が必要です*2。

カテゴリカルデータを変換する

カテゴリカルデータを数値として扱うために値をカラムとして変形します。

import pandas as pd

# csvファイルの読み込み
train_df = pd.read_csv("sanfranciscoCrime/input/train.csv")

# カテゴリカルデータを数値データに変換
dummyhours = pd.get_dummies(test_df.Dates.map(lambda x: pd.to_datetime(x).hour), prefix="hour", drop_first=True)
dummymonths = pd.get_dummies(test_df.Dates.map(lambda x: pd.to_datetime(x).month), prefix="month", drop_first=True)
dummyyears = pd.get_dummies(test_df.Dates.map(lambda x: pd.to_datetime(x).year), prefix="year", drop_first=True)
dummyweeks = pd.get_dummies(test_df.DayOfWeek, prefix="w", drop_first=True)
dummydistricts = pd.get_dummies(test_df["PdDistrict"], prefix="p", drop_first=True)

# train_dfを破壊的に更新する
train_df.drop(["Descript", "Resolution", "Dates", "DayOfWeek", "Address", "PdDistrict"], axis=1, inplace=True)
train_df = pd.concat([train_df, dummyyears, dummymonths, dummyweeks, dummyhours, dummydistricts], axis=1)

ランダムフォレストで予測

使用可能な分類器はランダムフォレスト以外にもサポートベクターマシンやロジスティック回帰、k近傍法などたくさんある。ランダムフォレストを選択した理由は、ある程度精度が高く、高速に分類できたためです。ランダムフォレストの分類器に与えたパラメータは、この後に述べるグリッドサーチで決めました。

from sklearn.ensemble import RandomForestClassifier

randomForest = RandomForestClassifier(n_estimators=100, max_depth=10, max_features=30)
y_train_df = train_df["Category"]
x_train_df = train_df.drop("Category", axis=1)

最適なパラメータを見つける

ランダムフォレストに与えるパラメータは、幾つか用意されています*3。

最適なパラメータを見つける手法の一つとして、グリッドサーチがあります。これを使うと、指定した全てのパラメータについて、モデルを生成し最適なパラメータを見つけてくれます。数学的な知識があればある程度どのようなパラメータを指定すればよいかわかるかもしれない。もしくは、時間と計算機資源があるならば地道にパラメータを探すのも良いかもしれない。ただし、今回のデータを用いてグリッドサーチを実行すると、組み合わせ総数の数だけモデルの生成とその精度の測定を実行しなければいけなくなるため、非常に時間がかかります。そのため、ここでは訓練データの1%のみしか用いていません*4。

from sklearn.model_selection import GridSearchCV

parameters = {
        'n_estimators'      : [20, 30, 50, 100],
        'max_features'      : [15, 20, 25, 30],
        'max_depth'         : [3, 5, 10, 15]
}
clf = GridSearchCV(RandomForestClassifier(), parameters)
clf.fit(x_train_df, y_train_df)

clf.grid_scores_
[mean: 0.19943, std: 0.00083, params: {'max_depth': 3, 'max_features': 5, 'n_estimators': 10},
 mean: 0.20057, std: 0.00191, params: {'max_depth': 3, 'max_features': 5, 'n_estimators': 20},
 mean: 0.20046, std: 0.00221, params: {'max_depth': 3, 'max_features': 5, 'n_estimators': 30},
 mean: 0.19920, std: 0.00064, params: {'max_depth': 3, 'max_features': 5, 'n_estimators': 50},
...
 mean: 0.23223, std: 0.00724, params: {'max_depth': 15, 'max_features': 20, 'n_estimators': 50},
 mean: 0.23850, std: 0.00771, params: {'max_depth': 15, 'max_features': 20, 'n_estimators': 100}]

その結果選ばれたパラメータは以下のとおりでした。場合によりますが、選ぶパラメータによって数％の精度の差が出るため、パラメータの選択は重要と言えるでしょう。

print(clf.best_params_)
print(clf.best_score_)
{'max_depth': 10, 'max_features': 30, 'n_estimators': 100}
0.242824601367

モデルの生成と分類

fit関数でモデルの学習を行い、predict関数でテストデータを分類します。ここでは、分類されたクラスを取得するpredictメソッドではなく、分類されるクラスの確率値を取得するpredict_probaメソッドを使用します。

randomForest.fit(x_train_df, y_train_df)
result = pd.DataFrame(randomForest.predict_proba(test_df), index=test_df.index, columns=randomForest.classes_)

結果

訓練データの規模でもモデルの学習までは行えましたが、レコード数があまりにも多かったためにテストデータの分類は実行できませんでした。おそらくディスク上に結果を書き出しつつ計算する方法もあるはずだと思いますが、今はそこまでは調べていません。ひとまず、小さなデータでは実行可能になったため、これで終わりとします。

その他のテクニック

このデータの分析を通して、学んだことを以降に記述します。

使用しない変数はdelメソッドで削除する

12MBのメモリを積んで学習と予測を実行してみたが、この規模のレコード数になるとOut of Memoryとなり、処理ができなかった。メモリの使用量を減らすことができるかどうかは、分析者にかかっている。delメソッドを使えば、明示的にメモリの解放もできる。

import gc

del train_df, dummyhours, dummyweeks, dummydistricts, dummymonths, dummyyears
gc.collect()

なるべく破壊的代入を行う

例えば訓練データからあるカラムを削除したいとするとき、新たにそのための変数を用意することはメモリを大きく消費することになります。具体的には、その変数に1GBの情報が詰められていたとするならば、ほんの少しだけ変更した変数を作成すると、前の変数の情報と合わせて2GBのメモリを消費することになります。dropメソッドにも"inplace"という破壊的な代入を行うための属性が存在するように、積極的に破壊的代入をしてべきでしょう*5。

train_df.drop(["Descript", "Resolution", "Dates", "DayOfWeek", "Address", "PdDistrict"], axis=1, inplace=True)

作成したモデルはディスクへ書き出す（うっかりデータが飛んでしまい、再計算する羽目になったから）

自分の動作環境があまり良くないためか、分析実行時にパソコンが制御不能に陥ることがよくありました。そのたびにモデルを学習し直さないといけなくなるととても厄介です。そういった時はsave&reloadができるようにしておいたほうがよいでしょう。

sklearn.externals import joblib
joblib.dump(randomForest, "sanfranciscoCrime/input/randomForest.pkl")
randomForest = joblib.load("sanfranciscoCrime/input/randomForest.pkl")

巨大なcsvファイルのデータを確認する

タイタニックのレコード数と比較すると、1,000倍規模になります。Excelを始めとするGUIツールを用いて巨大なcsvファイルを確認しようとすると、非常に時間がかかる。そのため、columnコマンドを使用すると、ファイルの中身を確認しやすい*6。

column -s, -t < train.csv | less -#2 -N -S

多重共線性を排除する

kaggleの幾つかのカーネルを見ていたとき、カテゴリカル変数の値を一つのカラムとして扱う際に、1つだけカラムを削除しているのをよく見かけました。これは多重共線性*7が発生することを回避するためであるようです*8。カテゴリカル変数からダミー変数を作り出すget_dummyメソッドにはdrop_first属性が与えられており、Trueのとき、最初の変数を削除してくれます。

dummyyears = pd.get_dummies(test_df.Dates.map(lambda x: pd.to_datetime(x).year), prefix="year", drop_first=True)

少量のデータで実験する

今回のデータは大きかったため、一度訓練データをtrain_dfに読み込んだ後に、サンプリングをしました。少量のデータで実験する場合にはサンプリングをするのもよいでしょう。

train_df = train_df.sample(frac=0.01)

*1:https://www.kaggle.com/maite828/test-random-forest

*2:もしかしたら、時間に関連するものは時系列データとして扱うことで処理ができるのかもしれない。

*3:3.2.4.3.1. sklearn.ensemble.RandomForestClassifier — scikit-learn 0.19.0 documentation

*4:そのため、精度が20%程度とかなり低い値になっている

*5:昔に関数型言語で実装していたため、副作用を発生させる文法を推奨するのには違和感がある

*6:エスケープ処理に対応していなかったためこれでも不十分だったが、純粋にlessコマンドを使用して開くよりはよい

*7:多重共線性とは？〜概要と対応方法〜 | 株式会社サイカ

*8:実際に、N個のダミー変数を作らずともN-1のダミー変数で全てを表現できる

2017-08-19

ベイズ主義的な解法による事後確率の更新

はじめに

ベイズの定理とは、ある変数について知られている確率（事前確率）と事象が起きたあとに得られる確率（事後確率）の関係を示す定理です。データサイエンスの業界で広く応用されています。

アメリカのクイズ番組を題材としたモンティ・ホール問題と呼ばれる確率の問題でちょっとした騒動が起きたそうですが、ベイズの定理を使えば簡単に解くことができます。しかし、ポール・エルデシュという数多くの論文を世に公表した高名な数学者であっても、その解法をうっかり間違えてしまったそうです*1。このように論争が起こる程度には、確率という分野は奥が深いようです。

この記事でモンティ・ホールの問題を解きませんが、代わりにクッキー問題の例を通してベイズの定理を紹介します。

頻度主義とベイズ主義

真の分布*2が存在することを仮定し、そこからデータが抽出される（事象が発生する）という考えに基づいた古典統計学（頻度主義）と、真の分布は存在せず事前知識と観測した事象によって構成される確率分布から確率が定まるという考えに基づくベイズ統計（ベイズ主義）の２つに分かれています。近年では、計算機の性能が向上したことによってビッグデータを解析する環境が整ったことや、膨大なパラメータ数を持つデータの収集が可能になったことから、サンプル数が少なくても良い精度を求めることができるベイズ統計がよく研究・応用されているようです。

この記事の目的

インターネット上では、ベイズの定理に関する記事は多く見受けられます。ただし、確率の概念といった抽象的なテーマについて論じるものが多かったです*3。

そして、事後確率がどのように更新されていくのかを具体的に述べたものが見つかりませんでした。より直感的に理解するために、ベイズの定理の更新に焦点を当てて説明します。

ベイズの定理とは

まずはベイズの定理の数理的性質を説明します。ベイズの定理*4は乗法定理とその対称性から導出されます。

$\displaystyle p\left( H | D\right) = \frac{p\left( D | H\right) p(H)}{p(D)}$

$H$ は仮説（hypothesis）を意味し、 $D$ はデータ（data）を意味しています。つまり、左項が示す意味は、"データが与えられたときの仮説の事後確率"となります。

クッキー問題とは

ベイズの更新を考えるために、クッキー問題を考えてみます。この問題は、"オライリー社のThink Bayes プログラマのためのベイズ統計入門"*5で紹介されていますが、より詳しく説明してみます。下図に示すように、2種類のクッキー（バニラ・チョコ）が複数枚入った2つのボウルがあります。

種類（枚）	ボウル1	ボウル2
バニラ	30	20
チョコ	10	20

ランダムにボウルを選択し、その中からランダムにクッキーを一つ取り出すことを考えます。この問題では、選ばれた種類がバニラであるとき、それがボウル1である確率を求めることを目的とします。求めたいのは、"ボウル1であるときにバニラを選ぶ確率"（ $\frac{30}{40}=0.75$ ）ではありません。なお、簡単のため、クッキーがボウルから減っていくことは考えないこととすることに注意してください。

クッキー問題をベイズの定理で解く

このとき、仮説 $H$ を"バニラをボウル1から選択する"、仮説 $D$ を"ボウル1からクッキーを選択する"*6とします。ベイズの定理における数式の意味は下表の通りです。

数式	役割	意味
$p\left( H \vert D\right)$	事後確率	バニラがボウル1から選ばれる確率
$p\left( D \vert H\right)$	尤度	ボウル1からバニラを選ぶ確率
$p(H)$	事前確率	ボウル1からクッキーを選ぶ確率
$p(D)$	正規化定数	２つのボウルからバニラを選ぶ確率

それぞれの役割を順番に説明します。

事後確率

事後確率 $p\left( H \vert D\right)$ は、データが与えられたときの仮説の確率を指します。ベイズの定理から事前確率と尤度、正規化定数の逆数の積で求められます。

尤度（関数）

尤度 $p\left( D \vert H\right)$ は、仮説が与えられたときのデータの確率を指します。例えば、ボウル1からクッキーを選ぶという行動を起こすとき、バニラというデータが得られる確率は、そのボウルの中のバニラの比率である0.75となります。

事前確率

事前確率 $p(H)$ は、仮説が起こる確率を指します。今回は2つのボウルをランダムに選択することから、それぞれの事前確率に $\frac{1}{2}$ が与えられます。しかし、事前確率の初期値は主観的に定められるものです。例えば、ボウル1がとても綺麗な器であり、ボウル2が汚く衛生的でない器であるとするならば、事前確率の初期値を当確率にすることはできないでしょう。あくまでも本ケースでは、ボウル1とボウル2の設定が平等であるから簡単に決めることができたということに注意してください。

正規化定数

正規化定数 $p(D)$ は、事後確率を1に調整するために使われる定数です。事前確率と尤度が定まるとき、正規化定数を求められます。バニラがどのボウルから選ばれるかということに焦点を絞っているため、正規化しなければ尤度と計算した値の総和が1にならず、それは事後確率でなくなってしまいます。そうなってはこの値の意味を解釈したい私達にとっては都合が悪いです。よって、ある条件下で独立な仮説を全て網羅する（排他的かつ網羅性がある）とき、それが占める確率を正規化定数として除算することで値を都合よく規格化します*7。今回は、クッキーの種類がバニラという条件の下で、ボウル1、ボウル2から取り出したという２つの仮説があるため、仮説の確率と条件の確率の積の総和を表のとおりに求めます。

ベイズの定理では、尤度関数と事前確率の初期値の与え方で事後分布がどのように変化していくかが定まります。そのため、その初期値の与え方や尤度関数をきちんと設計することは非常に重要です。

行動	事前確率	尤度	正規化定数の逆数	事後確率
ボウル1からバニラを選択	$0.50$	$0.75$	$\frac{1}{0.50\times0.75+0.50\times0.50}$	$0.60$
ボウル2からバニラを選択	$0.50$	$0.50$	$\frac{1}{0.50\times0.75+0.50\times0.50}$	$0.40$

新たな事象が起きると確率が変化する

ベイズの定理に従うと、ボウル1からバニラを選択する仮説に基づいたバニラを選択する確率 $p\left( H \vert D\right)$ は $0.50\times 0.75\times \frac{5}{8}=0.6$ に改訂されます。このとき、新たな事象（データ）を見たために、確率が変化したことがわかります。ベイズによる更新では、1回目で導出した事後確率の値を2回目の事前確率として扱います。

行動	事前確率	尤度	正規化定数の逆数	事後確率
ボウル1からバニラを選択	$0.60$	$0.75$	$\frac{1}{0.60\times0.75+0.40\times0.50}$	$0.69$
ボウル2からバニラを選択	$0.40$	$0.50$	$\frac{1}{0.60\times0.75+0.40\times0.50}$	$0.31$

1 回目の事後確率は 2 回目の事前確率として与えられます。そして同様に計算をすると、2 回目の事後確率が計算できます。一般化すると、n 回目の事後確率を n+1 回目の事前確率として扱うことで、ベイズによる事後確率の値は更新されていきます。下図はあらゆる分岐パターンを網羅した計3回の事象が発生したというデータを与えられたときに、事象発生地点で確率がどのように推移していくのかを示した図です。

f:id:GloSta:20170901020542p:plain

ボウル1からバニラを選択するという事象を引き当てると、徐々にボウル1からバニラを選ぶ確率 $p\left( H \vert D\right)$ が高くなり、ボウル1からチョコを引き当てると、その逆になることがわかります。ちなみに頻度主義に基づいて3回ともボウル1からバニラを選択した場合、その確率は100%という極端な値になってしまいます。そのため、データ数が明らかに少ないような本ケースの場合はベイズ統計のほうが好まれます。

まとめ

ベイズ統計を用いたときに確率がどのように推移していくのかを示しました。頻度主義的に問題を解くか、ベイズ主義的に問題を解くかは、状況しだいです。使い分けて、ベイズ統計の良い部分を活用していきましょう。

*1:

頻度確率による論理的な判断は一流の数学者でも間違えることがある？｜馬場正博の「ご隠居の視点」【寄稿】 | GiXo Ltd.

*2:大体の場合は真の分布の形はわからない

*3:個人的に最も理解しやすいベイズの定理の抽象的なお話はこちら

[教材] 今更だが, ベイズ統計とは何なのか. - ill-identified diary

*4:両辺に右辺の分母を掛け算すると左辺と右辺は事象 $H$ と事象 $D$ の同時確率になるので、覚えやすい

*5:https://www.amazon.co.jp/Think-Bayes-%E2%80%95%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%9E%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%85%A5%E9%96%80-Allen-Downey/dp/4873116945

*6:以降、この仮説をデータとして考えるとわかりやすい

*7:だからこそ、正規化定数を言葉として意味を与えるには難しいことが多い...?