はじパタ全力解説: 第10章 クラスタリング

[「はじめてのパターン認識」](https://amzn.to/2JlEWRz)の第10章「クラスタリング」の解説です。
前半(10.3まで)は、部分空間法と変わり難しい数式も出てこず理解しやすかったです。10.4 確率モデルによるクラスタリングで数式が出てきて難しくなり、ついに**EMアルゴリズムは理解を放棄しました**:cold_sweat: といっても、理解できなかったというより、時間かかるわりに、使うかの自信がなかったのでほとんど読まず飛ばしました。タイトルの「全力」に恥じる行為ですがお許しください。
※はじパタ勉強系は[記事「文系社会人がはじパタで機械学習を数式含めて理解した奮闘記」](https://qiita.com/FukuharaYohei/items/c56eaaa06bc85237f2f5)にまとめました。

# 参考リンク

- [はじパタ 10章 クラスタリング 前半(Slideshare)](https://www.slideshare.net/yamakatu/hajipata10)
- [はじめてのパターン認識輪読会 10章後半(Slideshare)](https://www.slideshare.net/kojiono507/10-31897824)
- [はじめてのパターン認識 10章](http://2gmon.github.io/haji_pata/10/index.html): 数式参考

# 内容
## 10.1 類似度と非類似度
### 10.1.1 距離の公理
**クラスタリング**は、データやクラスタ間の類似度、または非類似度に従って、似たものを集めてクラスタを作る方法です。画像中の色や明るさなどの類似した性質を持つ領域を他領域と区別する**セグメンテーション**や**ベクトル量子化**など広い分野で利用されています。
データやクラスタ間の類似度を測る基本的な尺度が**距離**です。ある尺度d(x,y)が距離であるために、次の**距離の公理**を満たす必要があります。

1. 非負性: d(x,y) ≧ 0
2. 反射律: d(x,y) = 0 となるのは x = y であるときのみに限る
3. 対称性: d(x,y) = d(y,x)
4. 三角不等式: d(x,z) ≦ d(x,y) + d(y,z)
 
### 10.1.2 ミンコフスキー距離
距離尺度には、よく使われる**ユークリッド距離**があります。それ以外にも**ミンコフスキー距離**と呼ばれる一般的な距離の定義から、様々な距離を派生させることができます。
d次元のデータがN個あり、i番目のベクトルを$x_i=(x_{i1},\ldots,x_{id})^T$とすると、2つのベクトル$x_i, x_j$間の距離は以下の式です。

```math
d(x_,x_j)=( \sum_{k = 1}^d |x_{ik} - x_{jk}|^a )^{1/b}
```

このaとbの値を変えて様々な距離を生み出します。
**aの値が大きければ特徴間の差に大きな重みが加えられます**。
**bの値が大きければ差分累乗和に対する重みが小さくなります**。
以下がaとbの値を変えて派生した距離です。
※一般的にミンコフスキー距離はa=bで定義されていることが多いようです。また、a!=bだと距離の公理を満たさなくなることもあり、派生させる場合には注意が必要です。詳しくはコメント欄を参照ください。

|No.|a|b|距離|式|
|:-:|:-:|:-:|:-:|:-:|
|1|1|1|ミンコフスキー距離<br>(マンハッタン距離/市街地距離)|$d(x_i, x_j) = \sum_{k = 1}^d &#124;x_{ik} - x_{jk}&#124;$|
|2|2|2|ユークリッド距離|$d(x_i, x_j) = (\sum_{k = 1}^d &#124;x_{ik} - x_{jk}&#124;^2)^{\frac{1}{2}}$|
|3|2|1|ユークリッド距離の2乗|$d(x_i, x_j) = \sum_{k = 1}^d &#124;x_{ik} - x_{jk}&#124;^2$|
|4|∞|∞|チェビシフ距離|$d(x_i, x_j) =  \lim_{a \rightarrow \infty}(\sum_{k = 1}^d &#124;x_{ik} - x_{jk}&#124;^a)^{\frac{1}{a}}$<br>$= \max_k &#124;x_{ik} - x_{jk}&#124;$|

ミンコフスキー距離の派生ではないですが、距離の公理を満たすものとして**キャンベラ尺度**があります。データ正規化の仕組みが内包されています。原点に近いほど分母が小さくなり、重みづけされます。

```math
d(x_i, x_j) = \sum_{k = 1}^d \frac{|x_{ik} - x_{jk}|}{|x_{ik}| + |x_{jk}|}
```

ベクトル間の角度($\cos \theta$)を用いた**コサイン類似度**があります。[記事「言語処理100本ノック-87:単語の類似度」](https://qiita.com/FukuharaYohei/items/a11597c9dcfbbd212557)でやりました。
※方向余弦は、コサイン類似度のこと。「方向余弦」はベクトルだと思っていたので、混乱しました。2本のベクトルがどれくらい同じ向きを向いているのかを表す指標で、1に近いほど2本のベクトルは同じ向きに近いです。取り得る値は-1から1までであり、負の値があるので距離の公理は満たしていません。

```math
d(x_i, x_j) = \frac{x_i x_j}{|x_i| |x_j|} 
= \frac{\sum_{k = 1}^d x_{ik} x_{jk}}{\sqrt{(\sum_{k = 1}^d x_{ik}^2) (\sum_{k = 1}^d x_{jk}^2) }}
```

距離の種類はこの辺が非常にわかりやすかったです。

- [クラスター分析の手法①（概要）](https://www.albert2005.co.jp/knowledge/data_mining/cluster/cluster_summary)
- [このデータセットにはどの距離を用いればよいの？](https://datachemeng.com/use_of_distance/)
 
## 10.2 非階層クラスタリング(K-平均法)
非改装クラスタリングの代表的手法である**K-平均法(K-means法)**です。d次元のN個のデータ$D=\\{x_1, \ldots, x_N\\}, x_i \in \mathbb R^d$をデータ間の類似性(距離)を尺度にあらかじめ定めたK個のクラスタに分類します。記事[「K-means 法を D3.js でビジュアライズしてみた」](http://tech.nitoyon.com/ja/blog/2013/11/07/k-means/)は見るだけで概要が理解できるほど、すばらしい記事です。
[別記事「言語処理100本ノック-97(scikit-learn使用):k-meansクラスタリング」](https://qiita.com/FukuharaYohei/items/072272c076a7cb5d7830)でも使いました。

- 各クラスタの代表ベクトルの集合: $M = \\{ \mu_i, \ldots, \mu_K \\}$
- k番目の代表ベクトルが支配するクラスタ: $M(\mu_k)$
- i番目のデータがk番目の代表ベクトルが支配するクラスタ$M(\mu_k)$への帰属有無を表す帰属変数: $q_{ik}$(下式参照)

```math
{q_{ik} = \left\{
\begin{array}{ll}
1 & \boldsymbol{x_i} \in M(\boldsymbol{\mu_k})\\
0 & otherwize
\end{array}
\right.}
```

K平均法の評価関数を以下に設定。$q_{ik}$が1の場合に代表ベクトルとの距離を最小化させます。

```math
J(q_{ik}, \mu_k)= \sum_{i=1}^{N}\sum_{k=1}^{K}q_{ik}||x_i-\mu_k||^2
```

$\mu_{ik}$を偏微分により最適化します。

```math
\frac{\partial J(q_{ik}, \mu_k)}{\partial \mu_k}=2\sum_{i=1}^Nq{ik}(x_i-\mu_k)=0\\
⇒\mu_k=\frac{\sum_{i=1}^Nq_{ik} x_i}{\sum_{i=1}^N q_{ik}
}
```

最後の形を見ると、各クラスタの代表ベクトルは、そのクラスタに帰属するデータの平均ベクトルになっていることがわかります。$\mu_{ik}$を最適化しましたが、$q_{ik}$と同時に最適化が難しいので、次の手順で逐次的に最適化を実施します。

#### 0. 初期化
 N個のデータをランダムにK個のクラスタに振り分け、それぞれのクラスタの平均ベクトルを求め、$\mu_k(k=1, \ldots, K)$とします。

#### 1. q最適化
$\mu_k$を固定し、帰属変数$q_{ik}$を以下に従って決定。

```math
q_{ik} = \left\{
\begin{array}{ll}
1 & k = argmin_j ||\boldsymbol{x_i}-\boldsymbol{\mu_j}||^2\\
0 & otherwize
\end{array}
\right.
```

K-平均法の**収束先は初期値に依存するため、最適解に近い解を得るためには何回か初期値を変えて実行する必要があります**。各代表ベクトルが**支配する領域は、代表ベクトル間のボロノイ境界で決まります**。
代表ベクトルは平均ベクトルなので、データベクトルと一般に一致をしません。代表ベクトルをデータベクトルに限定した方法を**K-メドイド法(K-medoids法)**と言います。代表ベクトルに平均ベクトルを使うため、K-平均法はデータ間の距離を定義する必要があります。一方で、代表ベクトルにデータベクトルを使うK-メドイド法はデータ間の非類似度行列(距離の行列でもOK)が定義できれば計算できます。
メドイドは、クラスタ内のデータ点のうち、その点以外のクラスタ内の他の点との比類似度(距離でもOK)の総和が最初うになる点です。以下を満たす$x_i$が代表ベクトル$\mu_k$の領域のメドイドです。

```math
i = \arg \min_{x_i \in M(\mu_k)}\sum_{y \in \{M(\mu_k) - x_i\}} d(x_j, y)
```

K-平均法が距離の2乗で誤差が評価されるのに対して、K-メドイド法は1乗なので、**外れ値の影響が少ない**です。

#### 2. μの最適化
$q_{ik}$を固定し、セントロイド$\mu_k$を決める(上で紹介した偏微分の式)。

#### 3. 繰り返し
「1. q最適化」と「2. μの最適化」を収束するまで繰り返します。

## 10.3 階層クラスタリング(融合法)
クラスタリングされていないデータから、類似度が高い順に融合して徐々に大きなクラスタを作り、最終的に1つのクラスタに統合していく方法を**階層的クラスタリング**または**融合法**と言います。データ・クラスタを融合する過程を**樹状図(dendrogram)**と呼ばれる形で表現できます。
以下が融合法の手順です。
1. n = Nとします
2. n × n の距離行列を作ります
3. 最も距離が近い2つのデータやクラスタをまとめて、一つのクラスタにします
4. n = n - 1とします
5. n > 1 であれば2へ戻ります。n = 1 であれば終了します。

類似度を測る尺度は、ミンコフスキー距離、キャンベラ尺度などを使います。**クラスタ間の類似度**は、以下のような様々なものがあります。

参考: [階層的クラスタリング入門の入門](https://www.slideshare.net/maskot3/ss-121402523)

### 10.3.1 単連結法
2つのクラスタ間で最も類似度の高いデータ間の距離をクラスタ間の距離にする方法を**単連結法(single linkage method)**と言います(**最短距離法**とも呼ぶようです)。
以下が単連結法によるクラスタリングの特徴です。

- クラスタに1つデータが追加されると、他クラスタとの距離は小さくなるか、または変化しない(最短距離なので当然)。
- 大きなクラスタができる傾向があります。
- あるクラスタから同じ距離に2つのクラスタがある場合、どちらを選んでも結果は同じ(途中経過は少し変わるが、2つが最短距離なので次の融合で2つ目のクラスタが対象となrimasu
- )。
- 近いデータが別なクラスタに属する連鎖効果が表れる場合があります(近いデータから融合していくので蛇のように長細くなっていくイメージ)

### 10.3.2 超距離
2つのデータ$x_i, x_j$が融合する直前のクラスタ間の距離を**超距離(ultrametric)**と言い、$\widetilde{d}(x_i,x_j)$で表します。
※「超距離」は、超距離法のようなものがあるのかと混乱しましたが、距離の種類・呼び方についての解説の章です。わかりにくい・・・

### 10.3.3 完全連結法
単連結法が最近傍距離を基準にしたのに対して、**完全連結法(complete linkage method)**は**最遠隣距離**をクラスタ間の距離にする方法です。以下が特徴です。

- クラスタに1つデータが追加されと、他クラスタとの距離は大きくなる(単連結法の逆)か、または変換しない。
- 大きなクラスタになりにくく、同じようなサイズのクラスタができる傾向。
- 連鎖効果は現れない(遠いデータを選んでいくから)

### 10.3.4 群平均法
**群平均法(group average method)**は2つのクラスタ内の全ての**データ対間の距離の平均で**クラスタ間の距離を決めます。クラスタ内データの平均ベクトル間ではないのが注意点です。

### 10.3.5 ウォード法
**ウォード法(Ward method)**は、クラスタAとBの距離を、それらを融合したときの**クラスタ内変動の増分**で定義し、距離の小さなクラスタから融合します。クラスタ内変動の増分は以下の式で表します。

```math
\begin{align}
d(A, B) & = \sum_{x \in A, B} d(x, \mu_{AB})^2 - (\sum_{x \in A} d(x, \mu_A)^2 + \sum_{x \in B} d(x, \mu_B)^2) \\
& = S_{AB} - (S_A + S_B) 
\end{align}
```

上記の式の補足です。つまり、平均との距離の差(各クラスタ内のばらつきを考慮)をもって、融合の基準とします。

|式|説明|
|:-:|:--|
|$d(x,y)$|ユークリッド距離|
|$\mu_{AB}$|クラスタAとBを融合したクラスタの平均ベクトル|
|$\mu_A$|クラスタAの平均ベクトル|
|$\mu_B$|クラスタBの平均ベクトル|
|S|**変動**、つまり平均からの距離(偏差)の2乗||

**ウォード法は階層法の中でもっとも精度が高い**らしいです。
これ以外に**重心法**や**メディアン法**があります。

## 10.4 確率モデルによるクラスタリング
K-平均法や融合法によるクラスタリングは、1つのデータは1つのクラスタにのみ分類されるため、**ハードクラスタリング**と呼ばれます。ハードクラスタリングに対して、確率モデルを使い所属クラスタを確率的に決める**ソフトクラスタリング**があります。多くの確率モデルは単峰性の確率分布しか表現できないので、**複数の確率モデルの重み付け線形和**で全体の確率分布をモデル化します。クラスタ数をK、k番目のクラスタの確率モデルを$p_k(x)$とし、全体の確率分布を$p(x) = \sum_{k = 1}^K \pi_k p_k(x)$で表します。
$\pi_k$がk番目の**確率モデルの重み(混合比)**です。このような確率モデルを**混合分布モデル**と言います。確率モデルに正規分布を使う場合、**混合正規分布モデル(混合ガウス分布モデル)**と呼ばれます。混合正規分布モデルを使う場合、**最初にクラスタ数Kを決める必要があります**。
※[「はじめてのパターン認識輪読会 10章後半」(Slideshare)](https://www.slideshare.net/kojiono507/10-31897824)が非常に参考になり、見ながら理解していきました。

### 10.4.1 混合正規分布モデル
k番目のクラスタを表すd次元正規分布関数を以下で評価します。要は多変量正規分布モデルです。参考に1変数の正規分布を2行目に記載。
※N は「正規分布」を表す英語 "normal distribution" の頭文字から取られている

```math
\begin{eqnarray}
\mathcal{N}(x|\mu_k, \Sigma_k) & = & 
\frac{1}{(2\pi)^{1/d} |\Sigma|^{1/2}} \exp (- \frac{1}{2} (x - \mu_k)^T \Sigma_k^{-1} (x - \mu_k) )\\
\mathcal{N}(x|\mu, \sigma^2) & = & 
\frac{1}{(2\pi)^{1/2} |\sigma|^{1/2}} \exp (- \frac{1}{2} (x - \mu)^2 |\sigma|^{1/2} )
\ \ *1変数正規分布\\
\end{eqnarray}
```

全体の分布は、線形和なので、以下で表します。

```math
p(x) = \sum_{k = 1}^K \pi_k \mathcal{N}(x|\mu_k, \Sigma), \;\; 
0 \leq \pi_k \leq 1, \;\; 
\sum_{k = 1}^K \pi_k = 1
```

$\pi_k$を**混合比**と呼び、混合正規分布モデルのパラメータは以下で表します。

```math
\pi = (\pi_1, \ldots , \pi_K), \;
\mu = (\mu_1, \ldots , \mu_K), \;
\Sigma = (\Sigma_1, \ldots , \Sigma)
```

$\mu_k$はd次元ベクトル、$\Sigma_k$はd×d対角行列で (d+1)d/2個の異なる要素を持つので、全体で K(混合比)+dK(平均ベクトル)+(d+1)dK/2(共分散行列の対角分)個のパラメータを求める必要があります。
※[混合ガウスモデル (Gaussian Mixture Model, GMM) ](https://datachemeng.com/gaussianmixturemodel/) が全般的に参考になります。

### 10.4.2 隠れ変数と事後確率 
混合正規分布モデルでは、データからK組の**混合比、平均ベクトル、共分散行列**を推定するために、1つのデータがどのクラスタに属するかを推定する必要があります。1つのデータがK個のクラスタのどこに属するかを表現するK次元変数を使います。これはOne-Hot Encoding表現で属するクラスタ$z_k$が1となります。

```math
z = (z_1, \ldots , z_K)^T, \;\; 
\sum_{k = 1}^K z_k = 1, \;\;
z = (0, \ldots , 0, 1, 0, \ldots , 0)^T
```

この変数は、変数xが所属する**隠れたクラスタ**を指定しているため、**隠れ変数**と呼ばれます。もともと提供されたデータになく、正解ラベルを自分で設定する変数なので、隠れ、と呼ぶイメージでしょうか。
変数xと隠れ変数zの**同時分布**はベイズの定理を元いて以下に分解できます。

```math
p(x, z) = p(z)p(x \mid z)
```

この中で$p(z_k=1)=\pi_k$なので、**隠れ変数の分布p(z)**は以下のように表現できます。
※以下の式の$z_k$部分は$z_k$乗という意味ではないです。

```math
p(z) = \prod_{k=1}^K \pi_k^{z_k}\\
\because \pi_k^{z_k=0}=1
```

上記は理解に時間をかけたので、少し詳しく(正しい理解か少し自信なし)。

|式|説明|
|:-:|:-:|
|$K=2$|クラスタ数が2|
|$\pi=\(\pi_1,\pi_2)=(0.2, 0.8)$|クラスタ1への所属確率が0.2<br>クラスタ2への所属確率が0.8|
|$z=(1,0)$|クラスタ1に属している場合の隠れ変数<br>最終的に事後確率で最尤推定する|
|$p(z_1=1)=\pi_1^{z_1} \times \pi_2^{z_2}=0.2 \times 1$|クラスタ1に所属する確率が0.2|

観測データの**隠れ変数による条件付き分布**は以下の式。

```math
p(x | z) = \prod_{k = 1}^K (\mathcal{N}(x | \mu_k, \Sigma_k))^{z_k}\\
\because p(x | z_k = 1) = \mathcal{N}(x | \mu_k, \Sigma_k)
```

p(x)は同時分布p(x, z)をすべてのzについての総和となるので以下の式。

```math
\begin{eqnarray}
p(x)& = & \sum_{k = 1}^K p(z)\\
& = & \sum_{k = 1}^K \prod_{k = 1}^K \pi_k^{z_k} \prod_{k = 1}^K (\mathcal{N}(x | \mu_k, \Sigma_k))^{z_k}\\
& = & \sum_{k = 1}^K \pi_k \mathcal{N}(x|\mu_k , \Sigma_k)\\
\end{eqnarray}
```

これで、ようやく求めたかった**隠れ変数の事後確率**$\gamma(z_k)$が計算できます。

```math
\begin{eqnarray}
\gamma(z_k) = p(z_k = 1 | x)& = & \frac{p(z_k = 1) p(x | z_k = 1)}{p(x)}\\
& = & \frac{p(z_k = 1) p(x | z_k = 1)}{\sum_{j = 1}^K p(z_j = 1)p(x | z_j = 1)}\\
& = & \frac{\pi_k \mathcal{N}(x | \mu_k, \Sigma_k)}{\sum_{j = 1}^K \pi_j \mathcal{N}(x | \mu_j, \Sigma_j)}\\
\end{eqnarray}
```

### 10.4.3 完全データの対数尤度 
観測データと隠れ変数を併せて**完全データ**とし、**完全データの尤度を最大にするパラメータを求めます**。ただし、最尤推定値を直接求めることはできません。

### 10.4.4 Q関数 
完全データの**対数尤度の隠れ変数に関する期待値をとると、隠れ変数の事後確率**になります。隠れ変数の期待値を事後確率で置き換えた関数を**Q関数**と呼びます。

### 10.4.5 EMアルゴリズム, 10.4.6 Mステップの導出, 10.4.7 EMアルゴリズムの性質
「全力」と言いながら、省略します。なぜなら、理解に時間がかかりそうな反面、今後使うのか？と感じる部分もあるため、飛ばします。
必要となったらこのあたりのページを見て勉強したいと思います。

- [EMアルゴリズム徹底解説](https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb)
- [EMアルゴリズム(slideshare)](https://www.slideshare.net/sotetsukoyamada/em-36315279)
- [数式を使わずイメージで理解するEMアルゴリズム(slideshare)](https://www.slideshare.net/yag_ays/em-algorithm-animation)

# 「はじパタ」勉強系記事リンク

|章|項目|時間(h)|難易度|学んだこと|
|:-:|:--|:-:|:-:|:--|
|1|[はじめに](https://qiita.com/FukuharaYohei/items/728b3fc15b5ddc9d9762)|8.9|:star:|特徴の型, 特徴空間, 次元の呪い|
|2|[識別規則と学習法の概要](https://qiita.com/FukuharaYohei/items/0bae9bbaef0099627751)|12|:star:|ホールドアウト法,交差確認法 ,一つ抜き法 ,ブートストラップ法 ,バイアス・分散トレードオフ, 過学習|
|3|[ベイズの識別規則](https://qiita.com/FukuharaYohei/items/1f05e6550ed73329febd)|14.8|:star:|ベイズ識別規則, ROC曲線|
|4|[確率モデルと識別関数](https://qiita.com/FukuharaYohei/items/10979bbd934410e82c42)|18|:star::star:|平均ベクトル, 共分散行列, 標準化, 無相関化, 白色化, 正規分布, 最尤推定|
|5|[k最近傍法(kNN法)](https://qiita.com/FukuharaYohei/items/0314f73cebb471deb515)|8|:star::star:|最近傍法, ボロノイ境界, kNN|
|6|[線形識別関数(前半)](https://qiita.com/FukuharaYohei/items/d1b4300055584d3f1727)|30.4|:star::star:|正規方程式|
|6|[線形識別関数(後半)](https://qiita.com/FukuharaYohei/items/f1933c0430704f86c60a)|前半に時間は含む|:star::star::star:|フィッシャーの線形判別関数, 判別分析法, ロジスティック回帰|
|7|[パーセプトロン型学習規則](https://qiita.com/FukuharaYohei/items/410655034fdef5204879)|13.5|:star::star:|多層パーセプトロン, 誤差逆伝播法, シグモイド関数|
|8|[サポートベクトルマシン](https://qiita.com/FukuharaYohei/items/5375d9bfe6bd471288ea)|14.7|:star::star:|カーネルトリック, ν-SVM|
|9|[部分空間法](https://qiita.com/FukuharaYohei/items/47624bac931b2668e888)|15.4|:star::star::star:|主成分分析, 特異値分解, CLAFIC法, カーネル主成分分析, カーネル部分空間法|
|10|[クラスタリング](https://qiita.com/FukuharaYohei/items/c58f6eea336ce5d9c00f)|8.4|:star::star:|距離の公理, ミンコフスキー距離, K-平均法, 融合法, 混合正規分布モデル|
|11|[識別機の組み合わせによる性能強化](https://qiita.com/FukuharaYohei/items/502d2d2b99ccd3980d75)|11.6|:star::star:|ノーフリーランチ定理, 決定木, バギング, アダブースト, ランダムフォレスト|