機械学習

ガウシアンナイーブベイズを自作する

みなさん機械学習してますか?

ぼくはちょっと今とあるプロジェクトで機械学習を用いたウェブアプリを自作してるのですが、その際にガウシアンナイーブベイズを自作したので、これはいい機会だということでガウシアンナイーブベイズの解説記事を書いちゃおうと思います。

n番煎じだということは承知ですが、現在ベイズのあたりを勉強している誰かの役に立てばいいなって感じです。

ちなみに数式を使ってちゃんと説明します(宣言)


準備

以下では数式を使って説明していきたいのですが、ノーテーション(記法)を統一するためにここで定義しておきます。

まず、$d$個のレコードを持つデータセットを$D = (\boldsymbol{x_1}, \boldsymbol{x_2}, ..., \boldsymbol{x_d})$とします。

各レコードは$n$個の特徴量を持っていて、$\boldsymbol{x_i} = (x_i^1, x_i^2, ..., x_i^n)^T$というような形をしています。

今回は機械学習を使ってデータセットからのラベルの予測をするわけですが、レコード$\boldsymbol{x_i}$に対する正解のラベルを$y_i$と表します。

さて、本稿ではscikit-learnのirisのデータセットを使い、アヤメという花のラベル予測をします↓

from sklearn import datasets

iris = Datasets.load_iris()

print(iris.data)
###        結果↓
### array([[5.1, 3.5, 1.4, 0.2],
### [4.9, 3. , 1.4, 0.2],
### [4.7, 3.2, 1.3, 0.2],
### [4.6, 3.1, 1.5, 0.2],
### ...
### [5.9, 3. , 5.1, 1.8]])

print(iris.target)
### 結果↓
###array([0, 0, 0, 0, 0, 0, ..., 2, 2, 2, 2, 2])

ちなみにirisのデータセットは150個のレコードが入ってます。軽量データセットです。

irisのデータセットでは特徴量は4個, 150個のレコードが格納されていることになります。


ナイーブベイズの説明

ナイーブベイズって何ぞやということですが、簡単に言ってしまえば条件付き確率モデルです。

数式的に言えば、レコード$\boldsymbol{x_i}$が与えられた時、それがあるラベル$y$である確率を表すモデルで、

P(y \ | \ \boldsymbol{x_i}) = P(y \ | \ x_i^1, x_i^2, ..., x_i^n)

というように表せます。

ここでベイズの定理を使えば

P(y \ | \ \boldsymbol{x_i}) = \frac{P(y) P(\boldsymbol{x_i} \ | \ y)}{P(\boldsymbol{x_i})} = \frac{P(y, \boldsymbol{x_i})}{P(\boldsymbol{x_i})}

として変形させることができます。

ここで結合確率(joint probability)$P(y, \boldsymbol{x_i})$について、チェインルールにより

\begin{align}

P(y, \boldsymbol{x_i}) &= P(x_i^1, x_i^2, ..., x_i^n, y) \\
&= P(x_i^1 \ | \ x_i^2, ..., x_i^n, y) P(x_i^2, ..., x_i^n, y) \\
&= P(x_i^1 \ | \ x_i^2, ..., x_i^n, y) P(x_i^2 \ | \ x_i^3 ..., x_i^n, y) P(x_i^3, ..., x_i^n, y) \\
&= \ ... \\
&= P(x_i^1 \ | \ x_i^2, ..., x_i^n, y) P(x_i^2 \ | \ x_i^3 ..., x_i^n, y) ... P(x_i^{n-1} \ | \ x_i^n, y) P(x_i^n \ | \ y) P(y)
\end{align}

と分解することができます。

さて、ここまでは普通のベイズですが、ナイーブベイズがなぜナイーブかというと、それぞれの特徴量が独立という仮定をするからです。

つまり

P(x_i^k \ | \ x_i^{k+1}, ..., x_i^n, y) \approx P(x_i^k \ | \ y)

とするわけです。

よってナイーブベイズモデルで結合確率(joint probability)$P(y, \boldsymbol{x_i})$は

\begin{align}

P(y, \boldsymbol{x_i}) &\approx P(x_i^1 \ | \ y) P(x_i^2 \ | \ y) ... P(x_i^n \ | \ y) P(y) \\
&= P(y) \prod_{k=1}^n P(x_i^k \ | \ y)
\end{align}

として表されるわけです。

以上を整理して、ナイーブベイズモデルは

\begin{align}

P(y, \boldsymbol{x_i}) &= \frac{P(y) \prod_{k=1}^n P(x_i^k \ | \ y)}{P(\boldsymbol{x_i})} \\
&= \frac{1}{Z} P(y) \prod_{k=1}^n P(x_i^k \ | \ y)
\end{align}

ここで$Z$は$P(\boldsymbol{x_i)}=\sum_k P(y_k) P(\boldsymbol{x_i} \ | \ y_k)$のことで、スケーリングのためのパラメータです。(確率なのでモデルの全出力の和が1になるようにスケーリングする)

実際のモデルの運用に際して、レコード$\boldsymbol{x_i}$に対するモデルの出力$\hat{y_i}$は

\hat{y_i} = \arg\max_{k} P(y) \prod_{k=1}^n P(x_i^k \ | \ y)

となります。

これがナイーブベイズモデルです。


ガウシアンナイーブベイズモデル

先ほどの解説でナイーブベイズモデルは解説しましたが、ではガウシアンなナイーブベイズモデルとはなんぞやというと、ラベル$y$に対しての特徴量$x_i^j$の尤度がガウス分布であると仮定したナイーブベイズモデルです。

つまり

P(x_i^j \ | \ y) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left[ - \frac{(x_i^j - \mu)^2}{2 \sigma^2} \right]

とするモデルです。

ここで平均値$\mu$と分散$\sigma^2$は最尤法によって決定されます。

具体的に、

\mu = \frac{1}{m} \sum_{i=1}^m x_i^k \\

\sigma^2 = \frac{1}{m} \sum_{i=1}^{m} (x_i^k - \mu)^2

となります。

これらは正規分布をそれぞれ$\mu$と$\sigma^2$で偏微分することによって得ることができます。

ここで注意したいのは、$\mu$と$\sigma$は正解ラベルが一緒のものだけのレコードの集合によって計算することです。(markdownでのいい感じの記法ができませんでした)


実装

さて、ガウシアンナイーブベイズを理論的に解説したので今度はこれを実装しようと思います。

実装はPython + Numpyでいきます。

まずはガウシアンを計算する関数を用意します。

数式は

gaussian(x, \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left[ - \frac{(x - \mu)^2}{2 \sigma^2} \right]

でしたね。

これを実装するとこんな感じです。

def gaussian(x, mean, variance):

return np.exp(-np.square(x - mean) / 2. / variance) / np.sqrt(2. * np.pi * variance)

これは簡単ですね。ナメクジでも実装できます。

次に最尤法によって$\mu$と$\sigma^2$の値を計算しましょう。

数式は

\mu = \frac{1}{m} \sum_{i=1}^m x_i^k \\

\sigma^2 = \frac{1}{m} \sum_{i=1}^{m} (x_i^k - \mu)^2

でした。

実装すると

def gaussian_param(input_data, target_data):

input_num = input_data.shape[0] # データの個数
feature_num = input_data.shape[1] # 特徴量の個数
target_set = np.unique(target_data) # ラベルのパターン
target_num = target_set.shape[0] # ラベルのパターン数

# 各特徴量、ラベル毎のパラメータの計算
mean = np.zeros([target_num, feature_num])
variance = np.zeros([target_num, feature_num])
for label in target_set:
mean[label] = np.mean(input_[target_==label], axis=0)
variance[label] = np.var(input_[target_==label], axis=0)

return mean, variance

学習データでのラベルの出現確率についても計算しておきます。

def target_frequency(input_data, target_data):

# 学習データに対するラベルの出現確率を計算
target_frequency = np.zeros([target_num])
for label in target_set:
self.target_frequency[label] = np.sum(target_==label) / input_num

return target_frequency

それでは最後に得られたパラメータから予測をするモデルを実装しましょう。

frequency = target_frequency(iris.data, iris.target)

def predict(input_data, mean, var):
input_num = input_data.shape[0] # 入力データの数
target_num = mean.shape[0] # ラベルのパターン数

answer = np.zeros([input_num])
for ith_data in range(input_num):
candidates = np.zeros([self.target_num]) # 各ラベル毎の確率
for label in range(self.target_num):
candidates[label] = frequency[label] * np.prod(gaussian(input_data[ith_data], mean[label], var[label]))
answer[ith_data] = np.argmax(candidates)

return answer

だいたいこんな感じです。

ちなみにこれをsklearnみたいな感じのクラス構成にするとこんな感じになります。

class NaiveBayes:

def gaussian(self, x, mean, variance):
return np.exp(-np.square(x - mean) / 2. / variance) / np.sqrt(2. * np.pi * variance)

def fit(self, input_, target_):
input_num = input_.shape[0]
feature_num = input_.shape[1]
target_set = np.unique(target_)
target_num = target_set.shape[0]
self.target_num = target_num

self.target_frequency = np.zeros([target_num])
for label in target_set:
self.target_frequency[label] = np.sum(target_==label) / input_num

self.mean = np.zeros([target_num, feature_num])
self.variance = np.zeros([target_num, feature_num])

for label in target_set:
self.mean[label] = np.mean(input_[target_==label], axis=0)
self.variance[label] = np.var(input_[target_==label], axis=0)

return self

def predict(self, input_):
input_num = input_.shape[0]
answer = np.zeros([input_num])

for ith_data in range(input_num):
candidates = np.zeros([self.target_num])
for label in range(self.target_num):
mean = self.mean[label]
var = self.variance[label]
candidates[label] = self.target_frequency[label] * \
np.prod(self.gaussian(input_[ith_data], mean, var))
answer[ith_data] = np.argmax(candidates)
return answer

さて、モデルを動かしてみましょう。

ちなみにsklearnのナイーブベイズのページではsklearnのirisのデータ150個に対して正答率96%だそうです。

Screen Shot 2019-05-04 at 18.27.21.png

こちらの自作モデルを実行してみると

if __name__ == "__main__":

model = NaiveBayes()

from sklearn import datasets
iris = datasets.load_iris()

model = model.fit(iris.data, iris.target)
print("Number of mislabeled points out of a total 150 points : " ,end="")
print(np.sum(model.predict(iris.data)-iris.target != 0))
### -> The number of wrong prediction is : 6

となり、無事にsklearnと同様の結果を得られました。めでたし。


まとめ

本稿ではガウシアンナイーブベイズの理論的解説とそれの実装を行いました。

結果としてsklearnと同じ程度の結果を得ることができ、いい感じっぽいです。

ちなみにPythonでクラス内メソッドで返り値をselfにするというのを今回初めてやったのですが、こういう書き方もあるのかってちょっと驚きましたね。

機械学習は楽しいので皆さんもぜひ頑張っていきましょう。

あ、よければぼくのツイッターをフォローしてください!@komi_edtr_1230

それではお疲れ様でした〜