１６．混合ガウスモデルを定式化しよう - keep learning blog（キープラーニングブログ）

――どんな条件であれ、私には確信がある。神は絶対にサイコロを振らない。――

混合ガウスモデルを定式化しよう

混合ガウスモデルを定式化しよう

f:id:yuki0718:20191101025533j:plain

研究関係の記事ばかり続いて恐縮ですが、前回の音声認識に引き続き、今回は機械学習に関する記事を書きます。少し前の記事では、流行りのディープニューラルネットワーク（DNN）に的を絞って機械学習の基本的な考え方をまとめました。

keep-learning.hatenablog.jp

機械学習というのはコンピュータを教育する方法の総称であって、そこには公文式やZ会のような様々な流派が存在するのでした。確かにDNNは機械学習の世界を席巻していますが、古典的な手法がすべてDNNに置き換わるわけではなく、未だにクラシックな機械学習が現役の分野も多くあります。

そこで、今回は偉大な物理学者アインシュタインがこよなく憎んだ「確率論」*1によって定式化される、古典的な機械学習モデルをご紹介したいと思います。具体的なモデルを考える前に、まずはどういった課題にそのモデルが有用なのかという点からご説明します。

教師なし学習とクラスタリング

第11回の記事でご紹介した教師あり学習（Supervised Learning）による分類問題（Classification）では、ワインの香りや味などの特徴量を数値化したN個の教師データと正解ラベルによってニューラルネットワークモデルを学習しました。そして、得られた学習済みモデルに対して評価用データを入力すると、ワインの予測ラベルをOne-Hotベクトルの形で出力し、適切な産地に分類することができました。

keep-learning.hatenablog.jp

それでは、最初からどのワインも産地が分からない場合、すなわち、正解付きの教師データが手元にない場合に、ワインを産地ごとにグループ分けすることは不可能なのでしょうか。正解不明のまま本番一発勝負でデータをグループ分けする手法はないのでしょうか。

一般に、このような課題をクラスタリング（Clustering）と呼びます。クラスタリングは、教師なし学習（Unsupervised Learning）の一種であり、教師あり学習の一種である分類問題（Classification）とは親戚のような関係にあります。当然ながら、正解の情報がない分、クラスタリングの難易度は分類問題よりも高いと言えます。

クラスタリングでよく用いられる機械学習モデルとしては、特徴量データの正規分布（Normal Distribution）を仮定した、混合ガウスモデル（GMM: Gaussian Mixture Model）が有名です。以下、このGMMを実装するうえで必要な理論的バックボーンを、確率論の観点からご紹介していきます。

多変量正規分布（Multivariate Normal Distribution）

私は理系の学生なので、大学時代は統計学の授業などもいくつか取っていました。統計学を学ぶ上で避けては通れないのが正規分布や条件付き確率です。まずは正規分布（ガウス分布）の身近な例として、試験の得点分布を考えましょう。

周知のとおり、日本では試験の結果を得点以外に「偏差値」という数値で知らせる慣習があります。偏差値は、受験者の得点分布が正規分布であると仮定した場合に*2、その得点が受験者の中でどの位置にあるのかを示す数値です。

ここで、試験の得点分布に正規分布を仮定するということは、平均（Mean） $\mu$ と分散（Variance） $\sigma$ という２つのパラメータを既知としたうえで、試験で $x$ 点を取る人数が、以下の正規分布関数（Normal Distribution）に比例すると認めることを意味します。

$\displaystyle Norm(\boldsymbol{x}_{i} \| \mu, \sigma) \stackrel{\mathrm{def}}{\equiv} \frac{1}{(2\pi)^{1/2}}\frac{1}{(\sigma)^{1/2}}exp\left[-\frac{(x-\mu)^{2}}{2\sigma} \right]$

分散 $\sigma$ の平方根を標準偏差（Standard Deviation）と呼び、偏差値という数値は自分の得点と平均 $\mu$ との差を標準偏差で割って10倍し、50を加えたものです（つまり平均と等しい得点は偏差値50になる）。ここまでは高校数学で扱われている範囲です。

それでは、試験が国語・数学・英語のように複数ある場合には、どのように正規分布を定義すればよいのでしょうか。まず、試験が3種類あるため、試験の得点は3次元のベクトル、

$\boldsymbol{x}=\begin{pmatrix}x_{1} \\ x_{2} \\ x_{3} \end{pmatrix}$

で表せばよいと分かります。そして、得点と合わせて平均と分散もそれぞれ以下のようなベクトルと行列にする必要があるでしょう。

$\boldsymbol{\mu}=\begin{pmatrix}\mu_{1} \\ \mu_{2} \\ \mu_{3} \end{pmatrix}$

$\boldsymbol{\Sigma}=\begin{pmatrix}\sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{pmatrix}$

ここで、3種類の試験の得点分布に正規分布を仮定するということは、上述した平均ベクトル（Mean） $\boldsymbol{\mu}$ と共分散行列（Covariance） $\boldsymbol{\Sigma}$ という２つのパラメータを既知としたうえで、3種類の試験で $\boldsymbol{x}$ 点を取る人数が、以下の多変量正規分布関数（Multivariate Normal Distribution）に比例すると認めることを意味します。

$\displaystyle Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}, \boldsymbol{\Sigma}) \stackrel{\mathrm{def}}{\equiv} \frac{1}{(2\pi)^{d/2}}\frac{1}{(det\Sigma)^{1/2}}exp\left[-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{T} \Sigma^{-1} (\boldsymbol{x}-\boldsymbol{\mu}) \right] \quad$ （dは次元数3）

平均と分散というパラメータが、正規分布の形を決めています。上で「パラメータは既知」と簡単に言いましたが、統計学をかじったことのある人であれば、全受験者の得点データから平均ベクトルや共分散行列を計算するための公式をご存知なのではないでしょうか*3。

実は、これは一種の機械学習であり、全受験者の得点データに基づいて正規分布モデルの学習を行っていると理解することもできます。

言い換えれば、試験の得点分布として多変量正規分布を仮定した場合に、全受験者の得点データから平均と分散というパラメータを計算することは機械学習における「学習」と同義です。

最適なパラメータが決まることで学習済みモデルは完成します。そして、モデルの最適なパラメータを計算するためには、以下に説明する確率論の考え方が必要になります。

条件付き確率（Conditional Probability）

ある事象Aが起こったという条件の下で事象Bが起こる確率のことを条件付き確率（Conditional Probability）と呼びます。例えば、大雨が降った（事象A）という条件の下で、電車遅延が起こる（事象B）確率、のようなものです。大雨が降ったからといって必ず電車遅延が起こるわけではないため、100%にはならないでしょう。この条件付き確率 $P(B \| A)$ は、以下のような確率の比として定義されます。

$\displaystyle P(B \| A) \stackrel{\mathrm{def}}{\equiv} \frac{P(A \cap B)}{P(A)}$

ここで、 $P(A \cap B)$ は同時確率（Joint Probability）と呼ばれ、単に $P(A, B)$ と省略表記する教科書がほとんどです。上の例でいうと、大雨が降って且つ電車遅延も起こる場合の確率を意味します。 $P(A)$ は電車遅延の有無に関わらず大雨が降る確率であり、以下のようにすべての $B$ （今の例では $B$ が起こる場合と $B$ が起こらない場合）で同時確率 $P(A \cap B)$ の和をとれば計算できます。

$\displaystyle P(A)=\sum_{B} P(A \cap B)$

このように、多変数関数である同時確率に対し、ある変数について和をとるような数学的操作を周辺化（Marginalization）と呼びます。周辺化は、Excelの表計算でいう小計のようなものです。

さて一方で、形式的には2つの事象の立場を逆転した条件付き確率、すなわち、 $P(A \| B)$ も定義することができます。上の例で言うと、電車が遅延したという条件の下で、大雨が降っている確率を意味します。電車が遅延したからといって、必ずしも大雨が降っているとは限らない（人身事故など）ため、100%にはならないでしょう。この条件付き確率も同様に定義するなら、

$\displaystyle P(A \| B)=\frac{P(B \cap A)}{P(B)}$

となります。

ところで、ここで一つ気が付くことがあります。この式に出てくる $P(B \cap A)$ は、電車が遅れて大雨も降る同時確率を意味しますが、同時確率は条件付き確率のように「Aが起こった後にBが起こる」といったA, Bの順序に関係するものではないため、 $P(B \cap A)$ は $P(A \cap B)$ と全く同じ確率 $P(A, B)$ です。したがって、どんな事象A, Bであろうと、条件付き確率の定義により以下の式が両立しなければなりません。

$\displaystyle P(A \| B)=\frac{P(A, B)}{P(B)}$

$\displaystyle P(B \| A)=\frac{P(A, B)}{P(A)}$

両式から $P(A, B)$ について整理して等号で結ぶと、

$P(A \| B)P(B)=P(B \| A)P(A)$

両辺を $P(B)$ で割り算して、

$\displaystyle P(A \| B)=\frac{P(B \| A)P(A)}{P(B)}$

この式のことをベイズの定理（Bayes' TheoremまたはBayes' Rule）と呼びます。定理といっても、導出過程からしてあまりにも簡単すぎる内容で、単に条件付き確率の定義式を言い換えたものだと理解すれば十分です*4。

最尤推定法（Maximum Likelihood Estimation）

条件付き確率の考え方がなぜ機械学習と関係するのか、我々初学者にはいまいちピンと来ません（少なくとも私はそうでした）。このミッシングリンクを補完するのが、条件付き確率で定義される尤度（ゆうど、Likelihood）という概念です。

一般に、 $i$ 番目に観測された特徴量データ $\boldsymbol{x}_{i}$ が手元にあるとき、データに対応する予測値をパーセプトロンや正規分布などの適当なモデルで仮定するのが機械学習の出発点でした。モデルというのは何かしらのパラメータ $\theta$ （正規分布なら $\mu$ や $\Sigma$ 、パーセプトロンなら重み係数 $w_{i-j}$ ）によって規定されます。

確率論の立場では、パラメータが $\theta$ という値をとる確率を、手元にあるデータ $\boldsymbol{x}_{i}$ を条件とする条件付き確率 $P(\theta \| \boldsymbol{x}_{i})$ で表現します。この条件付き確率を直接何かの関数でモデル化したものは識別モデル（Discrimination Model）と呼ばれています。

しかしながら、一般的に識別モデルの直接的な関数表現を見つけることは困難なので、ここでは逆転の発想をします。すなわち、手元にあるデータ $\boldsymbol{x}_{i}$ で規定される関数に基づいてパラメータ $\theta$ が決まるという考え方ではなく、モデルのパラメータ $\theta$ は先に分かっていて、そのパラメータ $\theta$ で規定される関数からデータ $\boldsymbol{x}_{i}$ の方が生成されるという考え方です。これを生成モデル（Generative Model）と呼びます。生成モデルでは、特定のパラメータ $\theta$ を条件とする条件付き確率 $P(\boldsymbol{x}_{i} \| \theta)$ を直接何かの関数でモデル化します。

仮にデータが合計 $N$ 個あり、データの生成（観測）が確率的に独立な事象であるとすれば、そのような状態が実現する確率は、

$\displaystyle L=\prod_{i=1}^{N} P(\boldsymbol{x}_{i} \| \theta)$

という条件付き確率の積となります。この確率は、手元にあるデータがパラメータ $\theta$ の生成モデルによって再現されることのもっともらしさを表し、尤度（ゆうど、Likelihood）と呼ばれています。確率論の（１つの）立場では、手元にあるデータに最も矛盾がないパラメータ $\theta$ は、尤度を最大化するように選ばれた $\theta_{ML}$ だと考えることができます。

$\displaystyle \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \theta_{ML}=\begin{equation} \argmax_{\theta} \prod_{i=1}^{N} P(\boldsymbol{x}_{i} \| \theta) \end{equation}$

このようなパラメータ $\theta$ の決定方法を最尤推定法（Maximum Likelihood Estimation）と呼びます。なお、積の形だと扱いにくいので、実用的には尤度の対数を最大化するように $\theta$ を選びます（対数関数は単調増加関数なので最大値を与える引数 $\theta_{ML}$ は同じ）。

$\displaystyle \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \theta_{ML}=\begin{equation} \argmax_{\theta} \log{L}=\argmax_{\theta} \sum_{i=1}^{N} \log{P(\boldsymbol{x}_{i} \| \theta)} \end{equation}$

これで今回のモデルを説明するための材料は揃いました。

混合ガウスモデル（Gaussian Mixture Model）

ある $N$ 個の特徴量データ $\boldsymbol{x}_{i} (i=1,2, \cdots, N)$ が手元にあるとき、これらを $C$ 種類のクラス（ラベル）に分類することを考えます。どのデータも正解クラスが分からない場合、この問題は教師なし学習であり、クラスタリング問題です。そして、クラスタリングでよく用いられるのが、正規分布（ガウス分布）の重み付き線形結合からなる生成モデル、いわゆる混合ガウスモデル（GMM：Gaussian Mixture Model）です。

潜在変数と生成モデルの定式化

GMMによるクラスタリングでは、 $k (=1,2,\cdots,C)$ 番目のクラスに属するデータは $k$ 番目の正規分布から生成される確率が最も高いと仮定します。このインデックス $k$ は、分類したいクラスの番号であるとともに生成モデルの番号を表し、現実に観測することはできない仮想的な変数であることから、潜在変数（Latent Variable）と呼ばれています。言葉で説明しても分かりにくいので、確率論の立場から具体的に定式化してみましょう。

まず、GMMは生成モデルなので、パラメータ $\theta$ は既知という前提をとります。その条件下において、 $k$ 番目の潜在変数が選択される確率 $P(k \| \theta)$ を重み係数（混合係数）と呼び、以下のように定義します。

$\displaystyle P(k \| \theta)=w_{k} \quad$ （ただし、確率なので $\displaystyle \sum_{k=1}^{C} w_{k}=1$ を満たす。）

さらに、 $k$ 番目のモデルによってデータ $\boldsymbol{x}_{i}$ が生成される確率 $P(\boldsymbol{x_{i}} \| k, \theta)$ を、多変量正規分布で定義します。

$\displaystyle P(\boldsymbol{x_{i}} \| k, \theta)=Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})$

これでGMMが定式化できました。何かを仮定するのはここまでです。後はこのモデルから演繹的にパラメータ $\theta$ の計算まで進めることができます。

尤度の計算

生成モデルが定義できたら、次は尤度（ゆうど、Likelihood）を計算します。上述したように、尤度とは $N$ 個のデータに関する条件付き確率 $P(\boldsymbol{x}_{i} \| \theta)$ の積で与えられる量でした。この式は潜在変数 $k$ を含まない形をしているので、これを計算するためには上述した周辺化（Marginalization）という数学的操作を行う必要があります。すなわち、「 $k$ 番目のモデルが選ばれる確率 $P(k \| \theta)$ 」と「 $k$ 番目のモデルからデータ $\boldsymbol{x}_{i}$ が生成される確率 $P(\boldsymbol{x_{i}} \| k, \theta)$ 」との積からなる同時確率 $P(\boldsymbol{x}_{i}, k \| \theta)$ を、すべての潜在変数 $k$ について足し合わせて $k$ に関する周辺化を行います。

$\displaystyle P(\boldsymbol{x}_{i} \| \theta)=\sum_{k=1}^{C} P(k \| \theta)P(\boldsymbol{x_{i}} \| k, \theta)$

後はすべてのデータ $\boldsymbol{x}_{i}$ について条件付き確率の積をとり、

$\displaystyle L=\prod_{i=1}^{N} P(\boldsymbol{x}_{i} \| \theta)= \prod_{i=1}^{N} \sum_{k=1}^{C} P(k \| \theta)P(\boldsymbol{x_{i}} \| k, \theta)$

これで尤度が計算できました。

最尤推定法による公式の導出

ここからは、最尤推定法によってモデルを規定するパラメータ $\theta$ の推定を行います。最尤推定法とは、尤度の対数を最大化するような引数 $\theta$ を求めるという考え方でした。すなわち、

$\displaystyle \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \theta_{ML}=\begin{equation} \argmax_{\theta} \sum_{i=1}^{N} \log{\sum_{k=1}^{C} P(k \| \theta)P(\boldsymbol{x_{i}} \| k, \theta)} \end{equation}$

これにGMMの定義式を代入して、

$\displaystyle \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \theta_{ML}=\begin{equation} \argmax_{\theta} \sum_{i=1}^{N} \log{\sum_{k=1}^{C} w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})} \end{equation}$

求めるパラメータ $\theta$ は、 $w_{k}, \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}$ の３つです。最大値をとるようなパラメータを求めるには、arg maxの右側にある目的関数を各パラメータで偏微分して=0とした方程式を解けばよいと分かります。ただし、 $w_{k}$ には $\sum w_{k}=1$ という拘束条件があったので、ラグランジュの未定乗数 $\lambda$ を導入した目的関数、

$\displaystyle \sum_{i=1}^{N} \log{\sum_{k=1}^{C} w_{k}\frac{1}{(2\pi)^{d/2}}\frac{1}{(det\Sigma_{k})^{1/2}}exp\left[-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu_{k}})^{T} \Sigma_{k}^{-1} (\boldsymbol{x}-\boldsymbol{\mu_{k}}) \right]}-\lambda \left( \sum_{k=1}^{C} w_{k}-1 \right)$

を偏微分して=0の方程式を解きます。この先の計算は文章量の関係で省略しますが、大学初年度レベルの線形代数学と偏微分の知識があれば、以下のように尤度を最大化するパラメータの公式が求められます。

$\displaystyle w_{k}=\frac{\sum_{i=1}^{N} \gamma_{i,k}}{N}$

$\displaystyle \boldsymbol{\mu}_{k}=\frac{\sum_{i=1}^{N} \gamma_{i,k}\boldsymbol{x}_{i} }{\sum_{i=1}^{N} \gamma_{i,k}}$

$\displaystyle \boldsymbol{\Sigma}_{k}=\frac{\sum_{i=1}^{N} \gamma_{i,k}(\boldsymbol{x}_{i}-\boldsymbol{\mu}_{k})(\boldsymbol{x}_{i}-\boldsymbol{\mu}_{k})^{T} }{\sum_{i=1}^{N} \gamma_{i,k}}$

ただし、上式において $\displaystyle \gamma_{i,k} \stackrel{\mathrm{def}}{\equiv} \frac{w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum_{k=1}^{C} w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}$ と定義した。

ここで定義された $\gamma_{i,k}$ は、公式を見やすくするために便宜上定義したものです。しかし、実はこれが確率論的に意味を持つことが天下り的に分かります。すなわち、条件付き確率 $P(k \| \boldsymbol{x}_{i}, \theta)$ をベイズの定理によって変形したもの、

$\displaystyle P(k \| \boldsymbol{x}_{i}, \theta)=\frac{P(\boldsymbol{x}_{i} \| k, \theta)P(k \| \theta)}{P(\boldsymbol{x}_{i} \| \theta)}=\frac{w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum_{k=1}^{C} w_{k}Norm(\boldsymbol{x}_{i} \| \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}$

これは $\gamma_{i,k}$ の定義と完全に一致します。つまり、 $\gamma_{i,k}$ は単なる便宜上のものではなく、 $i$ 番目のデータが観測された（手元にある）という条件下で $k$ 番目の生成モデルが選ばれる確率 $P(k \| \boldsymbol{x}_{i}, \theta)$ を表しています。このような確率論的解釈が可能なため、 $\gamma_{i,k}$ は $i$ 番目のデータに対する $k$ 番目の生成モデルの負担率（Responsibility）と呼ばれています。

上述した公式によってパラメータ $\theta_{ML}$ を計算して学習済みモデルを得ることができれば、この負担率によってデータ $\boldsymbol{x}$ がどのクラスに属するのか判定できます。具体的には、 $\theta_{ML}$ と $\boldsymbol{x}$ を各クラス $k$ の負担率の式に代入して計算すると、その値は $\boldsymbol{x}$ がクラス $k$ に属する確率になります。つまり、GMMによるクラスタリングでは、そのデータの負担率が最も大きいクラス $k$ に属すると判定できます。

以上、公式が導かれたところで、収まりがいいので今回の記事はここまでにします。次回は、最尤推定法によって得られたパラメータの公式を使って、具体的にプログラムを実装し、クラスタリングを行う予定です。というかもう作りました。

今回ご紹介したGMMは、私のような初学者が確率論の導入編で習うような古典的モデルですが、計算がシンプルな割に汎用性が高いため、様々な分野で広く用いられています*5。

稚文をお読みいただきありがとうございました。

*1:ボーア・アインシュタイン論争という有名な書簡のやり取りがあります。科学史の名場面なのでぜひGoogleで調べてみてください。

*2:もちろん、多くの場合に試験の得点分布は厳密な正規分布になどなりません。単純に便利なのでそう仮定しているのが実情です。

*3:例えば、平均ベクトルはすべての受験者の得点ベクトルを足し算して受験者の人数で割った値です。そんなことは直観的に自明ではあるものの、この公式は最尤推定法という確率論によって「導出」されるものです。

*4:ベイズの定理は、当たり前すぎて何がすごいのか分からないと思いますが、ある確率を別の確率同士の積で表現できるところがミソになっています。数学の世界には、何度掛け算しても元の関数の定数倍になるような性質を持つ関数の組み合わせ（共役関数）が存在するので、こういう性質の確率関数をモデルとして仮定し、ベイズの定理を連続適用して積を繋げていくと、計算が定数の掛け算になって驚くほど簡単になります。ベイズ確率論の応用については、本題から横道に逸れてしまうので今回は割愛します。

*5:例えば、音声認識の分野では2000年にD. A. Reynoldsらが提唱したGMM-UBMモデルという話者認識モデルが有名です（http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.117.338&rep=rep1&type=pdf）。私も自分の研究でディープラーニングとの比較例として実装しました。