3A

クラスタリング

講義の方針

  • データ解析アルゴリズムは、生物学の研究に不可欠
  • R や python のパッケージとして、すぐダウンロードして使うことができる。
  • しかし、原理を知らずに使うと、アルゴリズムの出す結果を盲信することになり危険。
  • 正しい理解：結果を疑えるようになること。

階層的クラスタリング

  • クラスタリング
    • データ点をクラスタ（互いに類似したグループ）に分割する。
    • 距離尺度：ユークリッド距離、コサイン距離など。
  • 階層的クラスタリング
    • 最初は各点がクラスタとなる。
    • 二つのクラスタを選んで合併していく。
      • Single Link：二つのクラスタをまたぐ最近点の距離 chainingにより、ゆるいクラスタができる（非推奨）
      • Complete Link：二つのクラスタをまたぐ最遠点の距離 最も固いクラスタができるが、細かく分け過ぎる傾向がある。
      • Average Link：二つのクラスタをまたぐ点対の平均距離 最もバランスがとれており、通常はこれを用いる。
    • 最後は一つのクラスタになる。
    • デンドログラム：合併の過程を樹形図の形で表したもの。

K-Means

分割的クラスタリングの一手法。$K$ クラスにデータ点を分割する。

以下を繰り返す。

  • 各データ点 $\mathbf{x}$ を最も近い中心点 $\mu_k$ に割り当てて、クラスタ $c_1,\ldots,c_K$ を作成する。
  • 割当を元に中心点を計算しなおす。 $$\mu_k = \sum_{j\in C_i} \mathbf{x}_j$$

実装

In [1]:
import numpy as np
import matplotlib.pyplot as plt
Data.
In [2]:
np.random.seed(0)
In [3]:
N=30; K=5
vmin=0; vmax=1
In [4]:
# Training Data.
data = np.random.uniform(vmin, vmax, (N,2))
x,y = data.T
In [5]:
plt.scatter(x,y)
plt.title("Training Data.")
plt.show()
In [6]:
# Background Color
Xs,Ys = np.meshgrid(np.linspace(vmin, vmax, 100), np.linspace(vmin, vmax, 100))
XYs = np.c_[Xs.ravel(),Ys.ravel()]
Training.

自作のモジュール（kerasy.ML.KMeans.KMeans）を使って学習させます。

In [7]:
from kerasy.ML.KMeans import KMeans
In [8]:
model = KMeans(K=K, random_state=0)
In [9]:
model.fit(data)
Visualization.
In [10]:
Nfig = len(model.history)
col = 3
row = Nfig//col+1 if Nfig%col != 0 else Nfig//col
In [11]:
# For Visualization.
plt_model = KMeans(K=K)
In [12]:
fig = plt.figure(figsize=(6*col,4*row))
for i,(idx,mu) in enumerate(model.history):
    # Plot Training Data with respective class color.
    ax = fig.add_subplot(row,col,i+1)
    ax.scatter(x,y,c=idx)
    # Plot Representative points(µk) of each class.
    mux,muy = mu.T
    ax.scatter(mux,muy,marker='x',c="red",s=200)
    # Background.
    plt_model.mu = mu
    Z, _ = plt_model.predict(XYs)
    Zs = Z.reshape(Xs.shape)
    ax.pcolor(Xs, Ys, Zs, alpha=0.2)
    ax.set_xlim(vmin,vmax)
    ax.set_ylim(vmin,vmax)
    ax.set_title(f"Iteration.{i:>02}")
plt.show()

混合ガウス分布

混合ガウス分布は、次のように混合係数 $\pi_k$ でガウス分布を線形重ね合わせした形でかける。

$$p(\mathbf{x}) = \sum_{k=1}^K\pi_k\mathcal{N}\left(\mathbf{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k\right)\qquad (9.7)$$

ここで、二値ベクトル（1-of-K符号化法）$\mathbf{z}$ の周辺分布は、$0\leq\pi_k\leq1$ と $\sum_{k=1}^K\pi_k = 1 $ を満たす $\pi_k$ を用いて、

$$p\left(\mathbf{z}\right) = \prod_{k=1}^K\pi_k^{z_k}\qquad (9.9)$$

と書ける。また、$\mathbf{z}$ が与えられた時の $\mathbf{x}$ の事後分布は当然ながら次のガウス分布である。

$$p\left(\mathbf{x}|\mathbf{z}\right) = \prod_{k=1}^K\mathcal{N}\left(\mathbf{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k\right)^{z_k}\qquad(9.11)$$

EMアルゴリズム

観測データ $\mathbf{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_N\}$ に対する対数尤度関数は、$(9.7)$ から以下のように書ける。（※明示的にパラメータを記載しています。）

$$\ln p\left(\mathbf{X}|\boldsymbol{\pi,\mu,\Sigma}\right) = \sum_{n=1}^N\ln\left\{\sum_{k=1}^K\pi_k\mathcal{N}\left(\mathbf{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k\right)\right\}\qquad (9.14)$$

Expectation step

M stepで必要となる $\gamma(z_{nk})$ を求める。

この値は、$\mathbf{x}$ が与えられた下での $\mathbf{z}$ の条件付き確率であり、混合要素 $k$ が $\mathbf{x}$ の観測を「説明する」度合いを表す負担率(responsibility)として解釈することもできる。

$$ \begin{aligned} \gamma(z_{k})\equiv p(z_{k}=1|\mathbf{x}) &=\frac{p(z_{k}=1)p(\mathbf{x}|z_{k}=1)}{\sum_{j=1}^{K}p(z_{j}=1)p(\mathbf{x}|z_{j}=1)}\\ &=\frac{\pi_{k}\mathcal{N}(\mathbf{x}|\mathbf{\mu_{\rm{k}}}, \mathbf{\Sigma_{\rm{k}}})}{\sum_{j=1}^{K}\pi_{j}\mathcal{N}(\mathbf{x}|\mathbf{\mu_{\rm{j}}}, \mathbf{\Sigma_{\rm{j}}})}\qquad (9.13) \end{aligned} $$

Maximization step

M stepでは、尤度関数を最大化するように各パラメータの値を更新する。

※ 導出過程は混合ガウス分布のEMアルゴリズムに載せています。

$$ \begin{aligned} \boldsymbol{\mu}_k &= \frac{\sum_{n=1}^N\gamma(z_{nk})\mathbf{x}_n}{\sum_{n=1}^N\gamma(z_{nk})} & (9.17)\\ \boldsymbol{\Sigma}_k &= \frac{\sum_{n=1}^N\gamma(z_{nk})\left(\mathbf{x}_n-\boldsymbol{\mu}_k\right)\left(\mathbf{x}_n-\boldsymbol{\mu}_k\right)^{\mathrm{T}}}{\sum_{n=1}^N\gamma(z_{nk})} & (9.19)\\ \pi_k &= \frac{\sum_{n=1}^N\gamma(z_{nk})}{N} & (9.22) \end{aligned} $$

実装

Data.
In [13]:
np.random.seed(0)
In [14]:
N=150; K=3
In [15]:
# Train Data.
data = np.concatenate([
    np.random.multivariate_normal(mean=[0, 0], cov=np.eye(2), size=int(N/3)),
    np.random.multivariate_normal(mean=[0, 5], cov=np.eye(2), size=int(N/3)),
    np.random.multivariate_normal(mean=[3, 2], cov=np.eye(2), size=int(N/3)),
])
x,y = data.T
In [16]:
plt.scatter(x,y)
plt.title("Training Data.")
plt.show()
In [17]:
# Background Color
xmin,ymin = np.min(data, axis=0)
xmax,ymax = np.max(data, axis=0)
Xs,Ys = np.meshgrid(np.linspace(xmin, xmax, 100), np.linspace(ymin, ymax, 100))
XYs = np.c_[Xs.ravel(),Ys.ravel()]
Training.

自作のモジュール（kerasy.ML.MixedDistribution.MixedGaussian）を使って学習させます。

In [18]:
from kerasy.ML.MixedDistribution import MixedGaussian
In [19]:
model = MixedGaussian(K=K, random_state=1234)
In [20]:
model.fit(data, span=20)
Visualization.
In [21]:
Nfig = len(model.history)
col = 3
row = Nfig//col+1 if Nfig%col != 0 else Nfig//col
In [22]:
# For Visualization.
plt_model = MixedGaussian(K=K)
plt_model.N, plt_model.M = data.shape
In [23]:
fig = plt.figure(figsize=(6*col,4*row))
for i,(epoch,idx,mu,S,pi) in enumerate(model.history):
    # Plot Training Data with respective class color.
    ax = fig.add_subplot(row,col,i+1)
    ax.scatter(x,y,c=idx)
    # Background.
    plt_model.mu = mu
    plt_model.S = S
    plt_model.pi = pi
    Z  = plt_model.predict(XYs)
    Zs = np.sum(Z, axis=1).reshape(Xs.shape) # Likelihood.
    ax.pcolor(Xs, Ys, Zs, alpha=0.2)
    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax)
    ax.set_title(f"{epoch:>02}epochs.")
plt.show()
In [ ]:
hidden
Table of Contents
Published
Last Updated
Category
生物データマイニング論
Tags
Contact
Other contents
  • Home
  • Blog
  • Front-End
  • Kerasy
  • Python-Charmers
  • Translation-Gummy