Shint’s blog

ピアノ調律からデータサイエンスまで

MENU

k近傍法

k近傍法

k近傍法(k-Nearest Neighbors: k-NN)とは、既知のグループが複数ある時に、どのグループに属するかわからない新たなデータが、どのグループに属するのかを推定する手法です。

f:id:Shint:20210215230723p:plain

上図はk近傍法のイメージです。赤と青で示したデータがあった際に、新しいサンプルがどちらのグループに属するかを考えます。データの空間において新しいサンプルから近いデータを3つ選ぶと \, (k=3) \, 、多数決により新しいサンプルは青のグループに分類されます。

k近傍法のハイパーパラメタは次の2つです。

  • 距離
    • 数学的に「距離の公理」を満たす定義は多数あります(ユークリッド距離、マハラノビス距離など)。
  • 新しいサンプルの近傍点の個数 \, k \,
    • 上の図で \, k=3 \, では青に分類されますが、 \, k=7 \, では赤に分類されるように結果が異なり得ます。

上の例は分類問題の例でしたが、回帰にも使えます。

  • 分類問題の場合:近傍k個のサンプルのうち最も多いクラスを予測クラスとする
  • 回帰問題の場合:予測値として近傍k個のサンプルの平均値をとる

k近傍法の特徴としては、

  • ノンパラメトリックモデル
  • オンライン学習が容易:パラメータの更新がないため、新しいデータを加えるだけで予測が可能です

サポートベクトルマシン

サポートベクトルマシン

サポートベクトルマシン(Support Vector Machineの考え方を以下に示します。

f:id:Shint:20210222154524p:plain

二つのクラスに対し境界線を与える際に、サポートベクトルマシンは境界線と隣接する点の距離(マージン)が最大となるように境界線を決定します。この隣接するデータ点のことをサポートベクトルと言います(あるデータ点 \, (x_0, x_1, \cdots, x_n) \,は一つのn次元ベクトルとしてみることができるため「ベクトル」です)。

線形分離可能なデータに関して:ハードマージン

下図のように線形分離可能な場合の定式化をハードマージンといい、 以下、2次元の例を用いてマージンはどのように最大化されるか詳しく見てゆきます(n次元でも大筋は同様です)。

赤のサポートベクトルで示される点 \, S_{+} \,の位置を表すベクトルを \, \boldsymbol{X}_{S+} = ( X_{1S+},  X_{2S+})^{t} \,、そこから境界線 \, l \,におろした垂線の足 \, H_{+} \,の位置を表すベクトル \, \boldsymbol{X}_{H+} = (X_{1H+}, X_{2H+})^{t} \,とします。

まず、境界線 \, l \,

 
\begin{eqnarray*}
   w_1x_{1}+w_2x_{2} + b =0 \qquad \cdots \quad (1)  \\
\end{eqnarray*}


とすると、 \, H_{+} \, に対しては、

 
\begin{eqnarray*}
   w_1X_{1H+}+w_2X_{2H+} + b =0 \qquad \cdots \quad (2)  \\
\end{eqnarray*}


という関係が成り立ちます。

次に、境界線 \, l \, に対する直行するベクトルを考えます。 オレンジには、境界線 \, l \, の係数を成分とするベクトル \, \boldsymbol{w} \, を示しました。一方、 境界線 \, l \, は、 x_{1} \,  \, w_{2} \, だけ進むと、  x_{2} \,  \, w_{1} \, だけ減少する関係を示したものですから(図中の水色)、 図を見比べると、 \, \boldsymbol{w} \, は境界線 \, l \, と直行するベクトルということがわかります(大学の数学的には、これは(1)の左辺のgradientをとることでも容易に得られます)。

今、 \, \vec{H_+S_+} \,  \, \boldsymbol{w} \, は平行であるため、

 
    \begin{pmatrix}
        X_{1S+} - X_{1H+} \\
        X_{2S+} - X_{2H+}
    \end{pmatrix}
   =
   t
    \begin{pmatrix}
        w_1\\
        w_2
    \end{pmatrix}
   \qquad \cdots \quad (3) \\

ただし、tはゼロでない定数です。 以下、 \, \vec{H_+S_+} \, の大きさを求めるために、 \, t \, を計算します。

式(3)の両辺で \, \boldsymbol{w} \, との内積をとり

 
   w_1(X_{1S+}-X_{1H+})+w_2(X_{2S+}-X_{2H+}) =t(w_1^2+w_2^2)  \\

展開して、

 
\begin{eqnarray*}
   w_1X_{1S+}+w_2X_{1S+} + b \\
   - (w_1X_{1H+}+w_2X_{2H+} + b)
\end{eqnarray*}
 =t(w_1^2+w_2^2)  \\

左辺のカッコ内は式 \, (2) \, より0となるため、

 \
   t = \dfrac{w_1X_{1S+}+w_2X_{2S+} + b}{w_1^2+w_2^2} 
     = \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{S+} + b }{|\boldsymbol{w}|^2} \\

よって、

 
   \vec{H_+S_+} = \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{S+} + b }{|\boldsymbol{w}|}
    \begin{pmatrix}
        w_1/|\boldsymbol{w}| \\
        w_2/|\boldsymbol{w}|
    \end{pmatrix} \\

 \, \vec{H_+S_+} \, の大きさがサポートベクトルマシンにおけるマージンで、それを \, M \, とします。 また、Class赤に分類されるサポートベクトル以外のデータ点と境界線 \, l \, の距離は \, M \, より大きくなければいけませんから、 \, x_{i'} \, をClass赤に分類されるデータ点として、

 
    \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i'} + b }{|\boldsymbol{w}|} \geq M \\

が成立します。等号はサポートベクトルの点において成立します。 この関係性を保ちながら \, M \, を最大化することが、今、得たい結果となります。 境界式 \, l \, を動かす際のパラメタは、式(1)で表されるように \, \{  w_1, w_2, b \} \, または \, \{ \boldsymbol{w}, b \} \, となりますので、

 \displaystyle
   \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i'} + b }{|\boldsymbol{w}|} \geq M,
   \quad \max_{\boldsymbol{w},b } M
   \qquad \cdots \quad (4) \\

f:id:Shint:20210223170847p:plain

一方、クラス青に分類されるサポートベクトル \, S_- \, から境界線 \, l \, におろした垂線の足 \, H_- (X_{1H-}, X_{2H-}) \, とで得られるベクトル \, \vec{H_{-}S_{-}} \, は、境界線 \, l \, と反平行になりますので、 符号が逆になります。

 \displaystyle
    - \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i''} + b }{|\boldsymbol{w}|} \geq M,
   \quad \max_{\boldsymbol{w},b } M
   \qquad \cdots \quad (5)  \\

ここに、 \, x_{i''} \, はClass青に分類されるデータ点です。

Class赤では正、Class青では負の符号を表す \, y_i \, を導入すると、式(4)(5)はまとめて次のように書くことができます。

 \displaystyle
   y_i \dfrac{ \boldsymbol{w} \cdot \boldsymbol{X}_{i} + b }{|\boldsymbol{w}|} \geq M,
   \quad  \max_{\boldsymbol{w},b } M
   \qquad \cdots \quad (6) \\

ここで、 \, x_i \, は任意のデータ点です。

境界線 \, l \, の両辺を定数倍しても結果は変わらないため、 計算しやすいように、以下のような変換を行います。

 \displaystyle
   \tilde {\boldsymbol{w}} = \dfrac{\boldsymbol{w}}{M | w | },
   \quad
   \tilde {b} = \dfrac{b}{M | w | } \\

この変換の式。および、変換後のベクトルの大きさ \, | \tilde{w} | = 1/M  \, より、

 \displaystyle
   y_i ( \tilde{ \boldsymbol{w} } \cdot \boldsymbol{X}_{i} + \tilde{b} ) \geq 1,
   \quad  \max_{ \tilde{ \boldsymbol{ w }}, b } \dfrac{1}{ | \tilde{w} | } \\

2つめの条件について、最大値をとることは逆数の最小値をとることと等価ですから、

 \displaystyle
   y_i ( \tilde{ \boldsymbol{w} } \cdot \boldsymbol{X}_{i} + \tilde{b} ) \geq 1,
   \quad  \min_{ \tilde{ \boldsymbol{ w }}, b } \dfrac{ | \tilde{w} |^2 }{2} \\

というようにして、計算をします(2乗にしたり、1/2を掛けたりするのは、 最適化計算をしやすくしているだけなので、気にしないでOKです)。

これらの関数を最適化してゆくアルゴリズムに受け渡すことで 学習が進みます。

線形分離不可能なデータに関して:ソフトマージン

これまでの例は線形分離可能な場合を扱ってきましたが、 そうでない場合(線形分離不可能な場合)の定式化をソフトマージンといいます。

f:id:Shint:20210301211219p:plain

図のように、境界線から片側のマージン分だけ離れたところにできる境界を超平面といいます(ここでの例は「直線」ですが、3次元に一般化すると「平面」、3次元以上では「超平面」と呼びます)。 今は、Class 赤側を正(Class 青側を負)としていました。

線形分離不可能な場合、各超平面を超えたデータ(分離できないデータ)に対してペナルティ(損失)を与えます。

損失の与え方は幾通りも考えられますが、一番自然なのは上の図のように「超平面を超えたデータに対して、超平面からの距離を損失としてあたえる」かと思います。

具体的に考えると、

  • あるデータ \, X_i \, が超平面を超えない場合(  \, y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b) \leq M \, )、ペナルティは0
  • あるデータ \, X_i \, が超平面を超える場合( \, y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b) \gt M \, )、ペナルティは超平面からの距離

f:id:Shint:20210301210907p:plain
図.Class 赤に対するペナルティの与え方の例

となるので、まとめて、

 \displaystyle
   \xi_i = \max \left\{ 0, \, M - \dfrac{y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b)}{||w||} \, \right\} \\

のように書けます。結局、ソフトマージンSVMの行うことは、「マージン \, M \, を最大にしつつペナルティ合計を最小化するような境界線 \, l \, を与える」ということになります。

「最大」と「最小」が混在していますので、少々めんどくさいですね。 これは、ハードマージンの時と同じように変数変換をすれば、最小化問題に集約できます。 つまり、変数変換後のペナルティ \, \tilde{\xi}_i \,

 \displaystyle
   \tilde{\xi}_i = \max \left\{ 0, \, 1 - y_i (\boldsymbol{w} \cdot \boldsymbol{X}_i + b) \, \right\} \\

として、

 \displaystyle
   y_i ( \tilde{ \boldsymbol{w} } \cdot \boldsymbol{X}_{i} + \tilde{b} ) \geq 1 -  \tilde{ \xi}_i, \quad   \min_{ \tilde{ \boldsymbol{ w }}, b }  \left\{ \dfrac{ | \tilde{w} |^2 }{2} +  C \sum_{i=1}^{n} \tilde{\xi}_i  \right\} \\

となります。ここであらたに登場した \, C (\leq 0)\, は、ハイパーパラメタで二つの量のバランスを調整します。

以上の議論に出てきたペナルティ項をスラック変数といいます。

ロジスティック回帰とソフトマックス回帰

ロジスティック回帰

ロジスティック回帰(Logistic Regression)は、回帰という名前がついていますが、分類の用途で使われます。

確率 \, p \, であるカテゴリに属すことの起こりやすさを表す量として、オッズ比  \,  p/(1-p) \, があります。 そして、この対数をとった関数をロジット関数   \, L(p) \, といいます。

 
    L(p) := \log \dfrac{p}{1-p} \\

 \, L(p) \, は定義域 \, 0 \leq p \leq 1 \, で値域 \, -\infty \,  \, \infty \, となる関数です。 ここで、ロジット関数の値は、特徴量 \, \{ x_1, x_2, \cdots, x_n \} \, の線形結合で与えられると仮定します。

まず、特徴量の線形関係を(線形回帰の時と同様に)

 
\begin{eqnarray*}
    z &=& w_1 x_1 + w_2 x_2 + \cdots + w_n x_n + b \\
       &=& \boldsymbol{w} \cdot \boldsymbol{x} + b \\
\end{eqnarray*}

としますロジット関数がこのような特徴量の1次関数 \, z \, で記載されているとすると、

 
   \log \dfrac{p}{1-p} = z \\

これを \, p \, について解くと、

 
   p = \frac{1}{1+e^{-z}} \\

と、シグモイド関数となります。

f:id:Shint:20210205102639p:plain

実装としては、特徴量 \, \{ x_1, x_2, \cdots, x_n \} \, の線形結合をシグモイド関数に代入すればよいわけですから、 以下のようになります。

import numpy as np

def input(x, w, b):
    z = np.dot(x,w) + b
    y = 1/(1+np.exp(-z))
    return y

ロジスティック回帰の学習とは線形回帰同様に、重み \, \boldsymbol{w} \, 、バイアス \, b \,を適切に調整してゆくという事になります。それでは、コスト関数はどう考えるかというと、ここで確率統計で学んだ最尤法を用いてゆきます (確率統計を参照してください)。

シグモイド関数の形からわかるように、 \, p \geq 0.5 \, ならば「あるカテゴリに属す」、 \, p \lt 0.5 \, ならば「あるカテゴリに属さない」という分類ができます。 つまり、二値分類となりますので、


  t_i =\begin{cases}
    1 & (p \geq  0.5) \\
    0 & (p \lt 0.5) \\
  \end{cases}

となる \, t_i (i=1,2,\cdots, n) \,を導入すると、i番目の \, t_i \,  \, 1 \, \, 0 \,になる確率は、  \, p^{t_i} (1-p)^{1-t_i} \,とまとめてかけるので、n個分類した時に実現する同時確率 \,  L(\boldsymbol{w},b)  \,

 \displaystyle
   L(\boldsymbol{w},b) = \prod_{i=1}^{n} p_i^{t_i} (1-p_i)^{1-t_i} \\

となります。この対数をとったものが対数尤度となります。ここでは、(最大問題を最小問題として扱うため)マイナスをかけた値をとり、

 \displaystyle
   \log L( \boldsymbol{w}, b ) = - \sum_{i=1}^{n} \{ t_i \log p_i + (1-t_i) \log (1-p_i) \} \\

とします。この形の関数を交差エントロピー誤差(Cross Entropy Loss)関数と言います。

# クロスエントロピー誤差
def closs_entropy(p, t):
    loss = -np.sum( (t * np.log(p) + (1 - t) * np.log(1 - p)), axis=0)
    return loss

この関数を最小にするような \, p \, が「もっともらしい」値となります。

ソフトマックス回帰(多項ロジスティック回帰)

ロジスティック回帰は2つのカテゴリへの分類でした。 3つ以上の分類ができるように、シグモイド関数ソフトマックス関数に差し替えたモデルを、 ソフトマックス回帰といいます。

 \displaystyle
    p_i := \dfrac{e^{z_i}}{\sum_{j=1}^{n} e^{z_j} }  \quad where \quad z_i =  \boldsymbol{w_i} \cdot \boldsymbol{x} +b_i \\
# ソフトマックス回帰
def softmax(x, w, b):
    z = np.dot(w, x) + b
    p = np.exp(z) / np.sum(np.exp(z))
    return p

線形回帰と正則化

線形回帰

単純な回帰のモデルとして線形回帰があります。予測する値 \, y \, が、入力するn個の特徴量 \, \{ x_1, x_2, \cdots, x_n \} \, の一次関数で与えられると考えます。

 
    f(x_1, x_2, \cdots, x_n) = w_1 x_1 + w_2 x_2 + \cdots + w_n x_n + b \\

とします。 \, \{  w_1, w_2, \cdots, w_n \} \, は特徴量の係数で重みといい、また、 \, b \,は定数項でバイアスといいます。

ここで、 \, \boldsymbol{w} = (w_1, w_2, \cdots, w_n) \, および \, \boldsymbol{x} = (x_1, x_2, \cdots, x_n) \, と定義して、

 
    f(\boldsymbol{x}) = w_1 x_1 + w_2 x_2 + \cdots + w_n x_n +b = \boldsymbol{w} \cdot \boldsymbol{x} + b\\

のように表すと便利です。



あるコンビニのアイスクリームの販売個数 \, y \, を予測するタスクを考える。 アイスクリームは暑い日になれば売れると考えれば「気温」や「湿度」が特徴量になると考えられる。 これらの変数を \, x_1, x_2 \, と表記する \, (n=2) \,

次に、販売個数 \, y \,  \, x_1, x_2 \, の一次関数であると仮定する。 その際、「気温が1℃上昇」することと「湿度が1%上昇」することの売り上げに対する影響度は異なる。 それゆえ、その影響度を個別に調整する重み  \, w_1, w_2  \, を導入して、
 
    y = w_1 x_1 + w_2 x_2 + b \\
とする。ただし、 \, b \, は定数項である。


さて、モデルが学習するというのは、この重み \,  \{ w_1, w_2, \cdots, w_n \} \, が手元にあるデータを再現するように調整するということです。予測値があるデータに近いかどうかを表す指標として平均二乗誤差(Root Means Square Error:RMSE)が一般的です。

 \displaystyle
\begin{eqnarray*}
   RMSE &=& \sqrt{\dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i) } \\
             &=& \sqrt{\dfrac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] } \\
\end{eqnarray*}

 \, \{ y_1, y_2, \cdots, y_m \} \, は実際のデータで、 \, m \, はデータの個数です。 要は、この差を最小にするように重み \{ w_1, w_2, \cdots, w_n \} を調整すればよいということになります。

なお、RMSEは標準偏差と似ていますが、RMSEは各データと推定値との差です(次のサイトが参考になります:二乗平均平方根などの統計量について



この3日間の気温、湿度、アイスクリームの売り上げデータがあり、 これらから明日のアイスクリームの販売個数を予想するタスクでは、
  • データ数: m=3
  • 実際の販売個数: y_1, \, y_2, \, y_3
  • 予測販売個数: \displaystyle \, \boldsymbol{w} \cdot \boldsymbol{x}_1 +b, \,  \boldsymbol{w} \cdot \boldsymbol{x}_2 + b,  \, \boldsymbol{w} \cdot \boldsymbol{x}_3  +b \,


RMSEを最小化する重み \,  \{ w_1, w_2, \cdots, w_n \} \, を考える際には、平方根を中身だけを考えればよいので、 平均二乗誤差(Means Square Error:MSE)を扱うことが機械学習では一般的です。このように、学習を進めるにあたり最小化したい関数をコスト関数といいます。

 \displaystyle
\begin{eqnarray*}
   MSE &=& \dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i)  \\
             &=& \dfrac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] \\
\end{eqnarray*}

正則化手法

リッジ回帰

過学習は、訓練データの量やノイズに対してモデルの方が複雑すぎるときに起こり、モデルに制約を加えることで過学習を抑制するのでした(正則化)。

リッジ回帰(Ridge Regression)は線形回帰の正則化で、コスト関数に \, \alpha/2 \sum_{i=1}^{n} w_i^{2} \, という正則化を加えます。リッジ回帰のコスト関数を \, J(\boldsymbol{w}) \, とすると、

 \displaystyle
\begin{eqnarray*}
   J(\boldsymbol{w}) &=& \dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i)  + \dfrac{\alpha}{2} \sum_{i=1}^{n} w_i^{2}\\
             &=& \dfrac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] + \dfrac{\alpha}{2} \sum_{i=1}^{n} w_i^{2} \\
\end{eqnarray*}

このようにすることで、学習モデルは訓練データを再現するだけでなく、重み \,  \{ w_1, w_2, \cdots, w_n \} \, できる限り小さく保つという束縛条件が一つ加わることになります。

ここで、 \, \alpha \, は人間が設定する値(ハイパーパラメータ)です。 \, \alpha =0 \, とすればモデルはただの線形回帰となり、 \, \alpha \, を大きな値にすれば、重み \,  \{ w_1, w_2, \cdots, w_n \} \, はすべてがほとんどゼロになり、モデルの予測値はほぼ定数となります。

数学的には、正則化項としてベクトルの大きさを表すノルムを加えていることになります。リッジ回帰の場合は \, L^{2} \, ノルムとなります。


参考
 \, n \, 次元のベクトル \, \boldsymbol{x} = (x_1, x_2, \cdots, x_n) \, および、 \, p ( 1\leq p \leq \infty ) \, に対して、
 \,  
    \sqrt[p]{|x_1|^p+|x_2|^p+\cdots+|x_n|^p}
 \, \boldsymbol{x} \,  \, L^p \, ノルムといい、 \,|| \boldsymbol{x}||_p \, と表す。


ラッソ回帰

ラッソ回帰(Least Absolute Shrinkage and Selection Operator Regressio: Lasso Regression)はリッジ回帰と同様に、コスト関数に正則化項を加えますが、 \, L^{2} \, ノルムではなく、 \, L^{1} \, ノルムを使います。

 \displaystyle
\begin{eqnarray*}
   J(\boldsymbol{w}) &=& \frac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i)  + \alpha \sum_{i=1}^{n} |w_i| \\
             &=& \frac{1}{m} \sum_{i=1}^{m} [ ( \boldsymbol{w} \cdot \boldsymbol{x}_i+b) - y_i ] + \alpha \sum_{i=1}^{n} |w_i| \\
\end{eqnarray*}

ラッソ回帰は重要性の低い特徴量の重みを0にする傾向があります。言い換えれば、自動的に必要な特徴量を選択し疎なモデルとなるということです。

しかし、ラッソ回帰は特徴量の数が多い時や特徴量の間に強い相関があるときには挙動が不規則になることがあることが注意点です。

Elastic Net

ラッソ回帰の弱点を克服しようとするのがErasitc Netです。考え方は、リッジ回帰とラッソ回帰の合いの子です。

 \displaystyle
\begin{eqnarray*}
   J(\boldsymbol{w}) &=& \dfrac{1}{m} \sum_{i=1}^{m} (f(\boldsymbol{x}_i) - y_i)  
                                         + r \alpha{2} \sum_{i=1}^{n} |w_i|
                                         + \dfrac{1-r}{2} \alpha \sum_{i=1}^{n} |w_i|^{2} \\
 \end{eqnarray*}

 \, r=0 \, でリッジ回帰、 \, r=1 \, でラッソ回帰と等しくなります。

リッジ回帰で試した際に意味がある特徴量が疑われる際には、ラッソ回帰の代わりにElastic Netを使用する方が良いとされています。

確率統計学の参考書

完全独習統計学入門 (小島寛之

バリバリ計算ができるようになるわけではないですが、統計の考え方を勉強するには最適です。

完全独習 ベイズ統計学入門(小島寛之

非常にわかりやすく、面白い本です。

確率統計キャンパス・ゼミ 改訂6(マセマ)

安定のマセマ。丁寧に式が導出されているので、持っていると何かと重宝します。

入門 統計学(栗原伸一)

多くのことが、まとまって書いてあります。なんだかんだで統計検定受験時にこれとマセマが役に立ちました。

入門 信頼性工学(福井泰好)

半導体の信頼性評価をする際に、まずこの本で必要なことは勉強しました。E資格というよりは、統計で信頼性がどのように担保されているのかなどを知りたい方へ。

最尤法

例えば、まず「1」となる確率が \, p \, 、「0」となる確率が \, 1-p \, のベルヌーイ分布があり、  \, ( x_1, x_2, x_3, x_4, x_5 ) = (1,1,1,1,0) \,  \, n=5 \, の確率変数があるとしましょう。

確率Lは、

 \displaystyle 
    L(p) = p^4(1-p)

であらわされます。

「現実となった確率変数は、確率が最大のものが実現した」という仮定を最尤原理(principle of maximum likelihood)といいます。

ここで試しに \, p = 0.1  \, の場合を計算してみると、

 \displaystyle 
  L(p=0.1) = 0.1^{4} \times 0.9 = 9 \times 10^{-5} \\

今度は \, p = 0.9  \, の場合を計算してみると、

 \displaystyle 
  L(p=0.9) = 0.9^{4} \times 0.1 = 0.0651 \\

となり、後者の方が確率 \, L(p=0.9) \, が大きくなります。

つまりこの例では、「  \, ( x_1, x_2, x_3, x_4, x_5 ) = (1,1,1,1,0) \, という事象があった時、 「1」となる確率 \, p \, は、 ( \, p = 0.1  \, または \, p = 0.9  \, なら) \, p = 0.9  \, の方が「もっともらしい」 と言えます。

ただし、ここまで読んでくださった皆さんでしたら、「 \, p = 4/5 = 0.8 \, だろう、おいおい」と、 お思いになるかもしれません。では、その計算の根拠を得るために、もう少しこの例を一般化してみましょう。

最尤原理では、確率 \, L(p) \,  \, p \, のとりうる値における「もっともらしさ」を表す関数とみなします。 このような「もっともらしさ」を尤度(likelihood)、その関数を尤度関数(likelihood function)といいます。

すなわち、最尤法(maximum likelihood method)とは、尤度関数を最大にする \, p \, を考えることとなります。

上の例で \, L(p) \, の最大値を求めればよいのですから、

 \displaystyle 
    \dfrac{\partial}{\partial p}  L(p) = p^3(4-5p) = 0\\

で、 \, p=0.8 \, が最大値となります。先ほどの「5回中4回が「1」だから \, p=4/5 \, 」という推測できることに一致しています。 今まで、このように事象から個別の確率を求めていた際には、暗に最尤原理を仮定していたということだったんですね。

さらに一般化します。

確率分布の変数を \, \theta \, として、n個の確率変数 \, \{ x_1, x_2, \cdots, x_n \} \, の同時確率分布は \, \theta \, の関数とみなし \, f(x, \theta) \, とすると、尤度関数は確率関数の積となり、

 \displaystyle 
    L(\theta) = f(x_1, \theta) \, f(x_2, \theta) \, \cdots \, f(x_n,\theta) = \prod_{i=1}^{n} f(x_n, \theta) \\

となります。ただし、正規分布が( \, \mu, \sigma^{2} \, )の関数であるように、 \, \theta \, は必ずしも1つではありません。

尤度関数の形として、積よりも対数をとって和の形にした方が扱いやすいため、

 \displaystyle 
    \log L(\theta) = \sum_{i=1}^{n} \log f(x_i, \theta) \\

を考えるのが普通です(対数尤度)。現実に起こりそうな \, \theta \, とは、対数尤度を最大にする \, \theta \, ですから、

 \displaystyle 
   \dfrac{\partial \log L(\theta)}{\partial \theta} = 0 \\

の解を考えることによって求めます。

確率分布の例

ベルヌーイ分布

何かが起こった時に、起こる結果が「A or B」のように2つに分類できる試行のことベルヌーイ試行と言います。 例えば、サイコロを振った時に「事象A:5以上の目」「事象B:4以下の目」のようなケースです。

Aとなる確率を \, p \, 、Bとなる確率を \, q \, とし、

 
    P(A)  = p \qquad P(B)  = q = 1 - p

と表すことにします。

このようなベルヌーイ試行をn回行い、ちょうどk回だけ事象Aとなる確率は次のように表されます。

f:id:Shint:20210125221514p:plain

①は「 \, n \, 回中 \, k \, 回選択する組み合わせ」の意味でコンビネーション、②は「確率 \, p \,  \, k \, 回」、③は「確率 \, q \, \, n-k \, 回」を意味します。 このように、確率変数 X=k に対して( \, q=1-p \, を用いて)

 
    P(X)  =  {}_n \mathrm{C} _k q^k (1-p)^{n-k} \quad (k=0,1,2, \cdots, n)

であらわされる分布をベルヌーイ分布または二項分布といいます。

ベルヌーイ分布では次の関係性が成り立ちます。
 (1) 期待値   E [ X ] = np
 (2) 分散  V [ X ] = npq

(1)の証明

コンビネーションを階乗であらわし、 \, n \,  \, p \, で括り、二項定理を使用します。 また \, p+q=1 \, の関係も使います。

 
\begin{eqnarray*}
    E [ X ] &=& \sum_{k=0}^n k \cdot {}_n \mathrm{C}_k \cdot p^k q^{n-k} \\
    &=&   \sum_{k=0}^n k \dfrac{n!}{k!(n-k)!} \times p^k q^{n-k} \\
    &=&  np \sum_{k=1}^n \dfrac{(n-1)!}{(k-1)!(n-k)!} \times p^{k-1} q^{n-k} \\
    &=&  np \sum_{k=0}^{n-1} \dfrac{(n-1)!}{k!(n-k-1)!} \times p^{k-1} q^{n-k-1} \\
    &=&  np \sum_{k=0}^{n-1} {}_{n-1} \mathrm{C}_k \cdot p^{k-1} q^{n-k-1} \\
    &=&  np \cdot (p+q)^{n-1} \\
    &=&  np
\end{eqnarray*}


(2)の証明

分散 \, V [ X ] \, の関係式より、

 
\begin{eqnarray*}
V [ X ] &=& E [ X^{2} ] - E [ X ]^{2}  \\
    &=& E [ X^{2} ]  - E [ X ] + E [ X ] - E [ X ]^{2}  \\
    &=& E [ X^{2} - X ] + E [ X ] - E [ X ]^{2}  \\
    &=& E [ X(X-1) ] + E [ X ] - E [ X ]^{2}  \qquad \cdots \quad (1)  \\
\end{eqnarray*}


となりますから、まず、右辺第1項 \, E [ X(X-1) ] \, を計算します。上の式と同じ方針で、 \, n(n-1)p^{2} \, で括り、 残りの部分には二項定理を適用します。

 
\begin{eqnarray*}
E [ X(X-1) ] &=& \sum_{k=0}^n k(k-1) \cdot {}_n \mathrm{C}_k \cdot p^k q^{n-k} \\
    &=& \sum_{k=0}^n k(k-1) \dfrac{n!}{k!(n-k)!} \times p^k q^{n-k} \\
    &=& n(n-1)p^2 \sum_{k=2}^n \dfrac{(n-2)!}{(k-2)!(n-k)!} \times p^{k-2} q^{n-k} \\
    &=& n(n-1)p^2 \sum_{k=0}^{n-2} \dfrac{(n-2)!}{k!(n-k-2)!} \times p^{k-2} q^{n-k-2} \\
    &=& n(n-1)p^2 \sum_{k=0}^{n-2} {}_{n-2} \mathrm{C}_k \cdot p^{k-2} q^{n-k-2} \\
    &=&  n(n-1)p^2 \cdot (p+q)^{n-2} \\
    &=&  n(n-1)p^2 \\
\end{eqnarray*}

この結果を(1)に代入して、

 
\begin{eqnarray*}
V [ X ] 
    &=& n(n-1)p^2 + np - (np)^2 \\
    &=& np(np-p+1-np) \\
    &=& np(1-p)
\end{eqnarray*}


E資格では \, n=1 \, のケースが問われるようです。

ベルヌーイ分布
 確率変数X( x = 0,1)について、
  (1) 確率分布  P(X) = p^x (1-p)^{1-x}
  (2) 期待値  E [ X ] = p
  (3) 分散  V [ X ] = p(1-p)


マルチヌーイ分布

ベルヌーイ分布が「事象A or 事象B」のように2つに分類できる問題に対して利用できる分布であったのに対し、 3つ以上に分類できる問題にも対応できるように拡張した分布がマルチヌーイ分布です。

いずれか1つのみ「 \, 1 \, 」の値をとり、それ以外「 \, 0 \, 」の値をとる \, n \, 個のパラメタ{ \, x_1, x_2, \cdots, x_n \, }と、 それに対応する確率{ \,  p_1, p_2, \cdots, p_n \, } ( \, \sum_{k=1}^n p_k = 1 \, )について

 
    P(X=x)  =  \prod_{k=1}^N p_k^{x_k}

となる確率分布です。

例えば、サイコロの目「3」のみが \, X=1 \, それ以外が \, 0 \, 、それぞれの目が出る確率は \, 1/6 \, とすれば、 Xは、マルチヌーイ分布に従います。


正規分布

最もよく使われる分布に正規分布(またはガウス分布があります。

 
    N(x; \mu, \sigma^2) = \dfrac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\dfrac{(x-\mu)^2}{2\sigma^2} \right)


 \, \mu \, は平均、 \, \sigma^{2} \, は分散( \, > 0 \, )、 \, \sigma \, 標準偏差を表します。

 \, \mu = 0 \,  \, \sigma^{2} = 1 \, である場合、標準正規分布といいます。

 
    N(x; 0, 1) = \dfrac{1}{\sqrt{2 \pi }} \exp \left( -\dfrac{x^2}{2} \right)


任意の正規分布は標準正規分布に変換することができます。確率変数 \, X \, 正規分布 \, N(\mu, \sigma^{2}) \, に従うとき、 上の2つの式(の特に指数)を見比べて、変数変換 \, Z=\dfrac{X-\mu}{\sigma^{2}} \, をすれば、 \, N(0, 1) \, となることがわかるかと思います。