scikit-learnのソースコードリーディング（ナイーブベイズ分類）

個人的にはプログラミングの勉強は写経が一番頭に入る気がする、ということで読んでいた。

気になったところ

データに正規分布を仮定したときのナイーブベイズ分類器について。平均を$\mu$、分散を$\sigma^2$としたときの正規分布は

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

これのlogをとると、

$$ \begin{eqnarray} \log p(x;\mu, \sigma^2) &=& \log \{\frac{1}{\sqrt{2\pi \sigma^2}} \{\exp{-\frac{(x-\mu)^2}{2\sigma^2}}\}\} \\\
&=& -\frac{1}{2}\log (2\pi \sigma^2) - \frac{(x-\mu)^2}{2\sigma^2} \end{eqnarray} $$

ナイーブベイズ分類器の対数尤度関数は、データがK次元ベクトルで表現されていて、それがN個あるとすると、 $$ \begin{eqnarray} \log L(X, Y; \mu, \sigma) &=& \log(\prod_{n=1}^N p(\mathbf{x}_n, y_n))\\\
&=& \log(\prod_{n=1}^N p(y_n)p(\mathbf{x}_n|y_n))\\\
&=& \sum_{n=1}^N \log p(y_n) + \sum_{n=1}^N \log p(\mathbf{x}_n|y_n)\\\
&=& \sum_{n=1}^N \log p(y_n) + \sum_{n=1}^N \sum_{k=1}^K\log p(x_{nk}|y_n)\\\
&=& \sum_{n=1}^N \log p(y_n) + \sum_{n=1}^N \sum_{k=1}^K \{-\frac{1}{2}\log (2\pi \sigma_{y_nk}^2) - \frac{(x_{nk}-\mu_{y_nk})^2}{2\sigma_{y_nk}^2}\} \end{eqnarray} $$

サンプル$\mathbf{x}$に対して出力される予測ラベル$\hat{y}$は

$$ \begin{eqnarray} \hat{y} &=& \mathop{\arg,\max}\limits_y \log p(\mathbf{x}, y)\\\
&=& \mathop{\arg,\max}\limits_y \log p(y)p(\mathbf{x}|y)\\\
&=& \mathop{\arg,\max}\limits_y \{\log p(y) + \sum_{k=1}^K \{-\frac{1}{2}\log (2\pi \sigma_{yk}^2) - \frac{(x_k-\mu_{yk})^2}{2\sigma_{yk}^2}\}\} \end{eqnarray} $$

対数尤度関数をnumpyに落とすと

"""
sigma.shape = (n_classes, n_features)
mu.shape = (n_classes, n_features)
"""

joint_log_likelihood = []
for i in range(np.size(classes)):
    # 事前分布の対数
    log_prior = np.log(class_piror[i])
    # log p(x|y)の対数の初項
    log_gauss1 = -0.5 * np.sum(np.log(2 * np.pi * sigma[i, :]))
    # log p(x|y)の対数の第二項
    log_gauss2 = -0.5 * np.sum((X - mu[i, :]) ** 2 / sigma[i, :])
    # クラスiの尤度のlogを取った値
    joint_log_likelihood.append(log_prior + log_gauss1 + log_gauss2)

となる。と思っていた。ところがscikit-learnのGaussianNBの該当箇所を見て見ると、

    def _joint_log_likelihood(self, X):
        X = array2d(X)
        joint_log_likelihood = []
        for i in range(np.size(self.classes_)):
            jointi = np.log(self.class_prior_[i])
            n_ij = - 0.5 * np.sum(np.log(np.pi * self.sigma_[i, :])) # np.piの前に2がない
            n_ij -= 0.5 * np.sum(((X - self.theta_[i, :]) ** 2) /
                                 (self.sigma_[i, :]), 1)
            joint_log_likelihood.append(jointi + n_ij)