階層的クラスタリングと、任意の数のクラスタへの分割

2016年に作った資料を公開します。もう既にいろいろ古くなってる可能性が高いです。

（追記：新しい記事は <a href="https://qiita.com/maskot1977/items/a35ac2fdc2c7448ee526">階層的クラスタリングとシルエット係数</a> をご覧ください。）

本実習では教師なし学習の一種である階層的クラスタリングを行ないます。
    * 階層的クラスタリング とは何か、知らない人は下記リンク参照↓
        * [階層的クラスタリングとは](http://image.slidesharecdn.com/140914-intel-l-141004040412-conversion-gate02/95/-12-638.jpg) 
        * [クラスタリング (クラスター分析)](http://www.kamishima.net/jp/clustering/)

### まずはサンプルデータの取得から


```python
# URL によるリソースへのアクセスを提供するライブラリをインポートする。
import urllib
# ウェブ上のリソースを指定する
url = 'https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/iris.txt'
# 指定したURLからリソースをダウンロードし、名前をつける。
urllib.urlretrieve(url, 'iris.txt')
```




    ('iris.txt', <httplib.HTTPMessage instance at 0x104115cf8>)




```python
import pandas as pd # データフレームワーク処理のライブラリをインポート
df = pd.read_csv("iris.txt", sep='\t', na_values=".") # データの読み込み
```

### 得られたデータを確認します。全体像も眺めます。


```python
df.head() #データの確認
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Unnamed: 0</th>
      <th>Sepal.Length</th>
      <th>Sepal.Width</th>
      <th>Petal.Length</th>
      <th>Petal.Width</th>
      <th>Species</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>5.1</td>
      <td>3.5</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>4.9</td>
      <td>3.0</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>4.7</td>
      <td>3.2</td>
      <td>1.3</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>4.6</td>
      <td>3.1</td>
      <td>1.5</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>5.0</td>
      <td>3.6</td>
      <td>1.4</td>
      <td>0.2</td>
      <td>setosa</td>
    </tr>
  </tbody>
</table>
</div>




```python
df.iloc[:, 1:5].head() #特徴量のデータは２列目〜５列目
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Sepal.Length</th>
      <th>Sepal.Width</th>
      <th>Petal.Length</th>
      <th>Petal.Width</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5.1</td>
      <td>3.5</td>
      <td>1.4</td>
      <td>0.2</td>
    </tr>
    <tr>
      <th>1</th>
      <td>4.9</td>
      <td>3.0</td>
      <td>1.4</td>
      <td>0.2</td>
    </tr>
    <tr>
      <th>2</th>
      <td>4.7</td>
      <td>3.2</td>
      <td>1.3</td>
      <td>0.2</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4.6</td>
      <td>3.1</td>
      <td>1.5</td>
      <td>0.2</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5.0</td>
      <td>3.6</td>
      <td>1.4</td>
      <td>0.2</td>
    </tr>
  </tbody>
</table>
</div>




```python
# 図やグラフを図示するためのライブラリをインポートする。
import matplotlib.pyplot as plt
%matplotlib inline
from pandas.tools import plotting # 高度なプロットを行うツールのインポート
plotting.scatter_matrix(df[df.columns[1:5]], figsize=(6,6), alpha=0.8, diagonal='kde')   #全体像を眺める
plt.show()
```


![output_6_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/d13878a0-7b41-e027-b17e-5ccd7a11d6ab.png)


### さていよいよ階層的クラスタリングの実行です。
* 階層的クラスタリングを実行するにあたって検討すべき事項は３つ。
    * 「特徴量の定義」
    * 「距離の定義」(metric)
    * 「リンケージ手法」(method)

ここでは、特徴量の定義は既に済んでいる（　df[df.columns[1:5]　）ものとし、距離の定義とリンケージ手法の選択を行って階層的クラスタリングを実行します。


```python
# metric は色々あるので、ケースバイケースでどれかひとつ好きなものを選ぶ。
# method も色々あるので、ケースバイケースでどれかひとつ好きなものを選ぶ。
from scipy.cluster.hierarchy import linkage, dendrogram
result1 = linkage(df.iloc[:, 1:5], 
                  metric = 'braycurtis', 
                  #metric = 'canberra', 
                  #metric = 'chebyshev', 
                  #metric = 'cityblock', 
                  #metric = 'correlation', 
                  #metric = 'cosine', 
                  #metric = 'euclidean', 
                  #metric = 'hamming', 
                  #metric = 'jaccard', 
                  #method= 'single')
                  method = 'average')
                  #method= 'complete')
                  #method='weighted')
dendrogram(result1)
plt.title("Dedrogram")
plt.ylabel("Threshold")
plt.show()
```


![output_8_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/675c5cde-ad3b-f25c-8fca-5c8060ff18ca.png)


とりあえずこれで、デンドログラム（樹形図）は描けました。

### クラスタを得る

Threshold （閾値）を変えることでクラスタの数は変わります。どのくらいの閾値が適切か、どのくらいのクラスタ数が適切か、データの性質を理解した上でケースバイケースで検討しなければいけません。適切な閾値、あるいは適切なクラスタ数でクラスタを得る方法を作りましょう。


```python
# Threshold を変えるとクラスタ数や平均クラスタサイズがどう変わるか調べる
import numpy as np
import matplotlib.pyplot as plt

n_clusters = len(df)
n_samples = len(df)
df1 = pd.DataFrame(result1)
x1 = []
y1 = []
x2 = []
y2 = []
for i in df1.index:
    n1 = int(df1.ix[i][0])
    n2 = int(df1.ix[i][1])
    val = df1.ix[i][2]
    n_clusters -= 1
    x1.append(val)
    x2.append(val)
    y1.append(n_clusters)
    y2.append(float(n_samples) / float(n_clusters))
        
plt.subplot(2, 1, 1)
plt.plot(x1, y1, 'yo-')
plt.title('Threshold dependency of hierarchical clustering')
plt.ylabel('Num of clusters')
plt.subplot(2, 1, 2)
plt.plot(x2, y2, 'ro-')
plt.xlabel('Threshold')
plt.ylabel('Ave cluster size')
plt.show()
```


![output_10_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/c3fda50c-2b7c-7d0d-9344-6cf9bfa2f909.png)



```python
# 指定した thoreshold でクラスタを得る関数を作る
def get_cluster_by_threshold(result, threshold):
    output_clusters = []
    output_cluster_ids = []
    x_result, y_result = result.shape
    n_clusters = x_result + 1
    cluster_id = x_result + 1
    father_of = {}
    df1 = pd.DataFrame(result)
    x1 = []
    y1 = []
    x2 = []
    y2 = []
    for i in df1.index:
        n1 = int(df1.ix[i][0])
        n2 = int(df1.ix[i][1])
        val = df1.ix[i][2]
        n_clusters -= 1
        if val < threshold:
            father_of[n1] = cluster_id
            father_of[n2] = cluster_id
            
        cluster_id += 1
    
    cluster_dict = {}
    for n in range(x_result + 1):
        if not father_of.has_key(n):
            output_clusters.append([n])
            continue
        
        n2 = n
        m = False
        while father_of.has_key(n2):
            m = father_of[n2]
            #print [n2, m]
            n2 = m
        
        if not cluster_dict.has_key(m):
            cluster_dict.update({m:[]})
        cluster_dict[m].append(n)
    
    output_clusters += cluster_dict.values()
    
    output_cluster_id = 0
    output_cluster_ids = [0] * (x_result + 1)
    for cluster in sorted(output_clusters):
        for i in cluster:
            output_cluster_ids[i] = output_cluster_id
        output_cluster_id += 1
            
    return output_cluster_ids
```


```python
# 指定した thoreshold でクラスタを得る関数を使う。
print get_cluster_by_threshold(result1, 0.05)
```

    [0, 1, 1, 1, 0, 2, 1, 0, 1, 1, 0, 0, 1, 3, 2, 2, 2, 0, 2, 0, 0, 0, 4, 0, 0, 1, 0, 0, 0, 1, 1, 0, 2, 2, 1, 1, 0, 0, 1, 0, 0, 5, 1, 0, 2, 1, 0, 1, 0, 0, 6, 6, 6, 7, 6, 7, 6, 8, 6, 7, 8, 7, 9, 6, 7, 6, 7, 7, 9, 7, 10, 6, 10, 6, 6, 6, 6, 11, 6, 7, 7, 7, 7, 10, 7, 6, 6, 9, 7, 7, 7, 6, 7, 8, 7, 7, 7, 6, 8, 7, 11, 10, 12, 11, 11, 13, 14, 12, 12, 12, 11, 11, 11, 10, 10, 11, 11, 13, 13, 9, 11, 10, 13, 10, 11, 12, 10, 10, 11, 12, 12, 13, 11, 10, 10, 13, 11, 11, 10, 11, 11, 11, 10, 11, 11, 11, 10, 11, 11, 10]



```python
# 指定したクラスタ数でクラスタを得る関数を作る。
def get_cluster_by_number(result, number):
    output_clusters = []
    x_result, y_result = result.shape
    n_clusters = x_result + 1
    cluster_id = x_result + 1
    father_of = {}
    df1 = pd.DataFrame(result)
    x1 = []
    y1 = []
    x2 = []
    y2 = []
    for i in df1.index:
        n1 = int(df1.ix[i][0])
        n2 = int(df1.ix[i][1])
        val = df1.ix[i][2]
        n_clusters -= 1
        if n_clusters >= number:
            father_of[n1] = cluster_id
            father_of[n2] = cluster_id
            
        cluster_id += 1
    
    cluster_dict = {}
    for n in range(x_result + 1):
        if not father_of.has_key(n):
            output_clusters.append([n])
            continue
        
        n2 = n
        m = False
        while father_of.has_key(n2):
            m = father_of[n2]
            #print [n2, m]
            n2 = m
        
        if not cluster_dict.has_key(m):
            cluster_dict.update({m:[]})
        cluster_dict[m].append(n)
    
    output_clusters += cluster_dict.values()
            
    output_cluster_id = 0
    output_cluster_ids = [0] * (x_result + 1)
    for cluster in sorted(output_clusters):
        for i in cluster:
            output_cluster_ids[i] = output_cluster_id
        output_cluster_id += 1
            
    return output_cluster_ids
```


```python
# 指定したクラスタ数でクラスタを得る関数を使う。
print get_cluster_by_number(result1, 5)
```

    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 3, 2, 4, 2, 4, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 2, 2, 4, 4, 4, 4, 2, 4, 2, 4, 2, 4, 4, 2, 2, 4, 4, 4, 4, 4, 2, 2, 4, 4, 4, 2, 4, 4, 4, 2, 4, 4, 4, 2, 4, 4, 2]


以上で、階層的クラスタリングができるようになりました。でもこれだけでは、しっくり来る結果が得られたのかどうかよく分かりませんので、図示してみましょう。

### クラスタリング結果のマッピング
* クラスタリング結果を Scatter Matrix 上にマッピングします。
* 色は自分で決めることもできます（[総合実験４日目](http://nbviewer.jupyter.org/github/maskot1977/ipython_notebook/blob/master/%E7%B7%8F%E5%90%88%E5%AE%9F%E9%A8%93%EF%BC%94%E6%97%A5%E7%9B%AE.ipynb) 参照）が、色の数が増えると面倒なので次のように自動で決めることもできます。 [cmapのオプション](http://fififactory.com/wiki/index.php?plugin=attach&refer=FrontPage%2FPython%2FMatplotlib&openfile=mpl_cmap.png) を設定すればいろんな色の組み合わせで表現できます。


```python
# Scatter Matrix に、クラスタリング結果をマッピング
import matplotlib.pyplot as plt
import matplotlib.cm as cm
%matplotlib inline
from pandas.tools import plotting # 高度なプロットを行うツールのインポート
plotting.scatter_matrix(df[df.columns[1:5]], figsize=(6,6), alpha=0.8, diagonal='kde', 
                        c=get_cluster_by_number(result1, 5), cmap=cm.rainbow, color=[], s=50)
plt.show()
```


![output_17_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/7b8dd7ce-e086-23c2-9e43-6d97410c8d16.png)


* クラスタリング結果を、主成分分析のプロット上にマッピングします。


```python
# 主成分分析のプロットに、クラスタリング結果をマッピング
from sklearn.decomposition import PCA
pca = PCA()
pca.fit(df[df.columns[1:5]])
feature = pca.transform(df[df.columns[1:5]])
plt.figure(figsize=(6, 6))
plt.scatter(feature[:, 0], feature[:, 1], c=get_cluster_by_number(result1, 5), s=60, alpha=0.5,
            cmap=plt.cm.rainbow)
plt.show()
```


![output_19_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/7f29eeae-9a18-eba9-21e7-d4744e83c252.png)


### もうちょっと詳しく、流れを理解したい
以上で、階層的クラスタリングと、クラスタリング結果のマッピングまでできるようになりました。ここから先は、やらなくても良いと言えばやらなくても良いところです。階層的クラスタリングをするときに、何が起こっているのか少しだけ詳しく見ながら進めると次のようになります。

__距離行列生成__


```python
# 距離といっても色んな距離があるので、どれか選ぶ。
import scipy.spatial.distance as distance
n_samples = len(df)
dMatrix = np.zeros([n_samples, n_samples])
for i in range(n_samples):
    for j in range(n_samples):
        dMatrix[i, j] = distance.braycurtis(df.ix[i, 1:5], df.ix[j, 1:5]) # braycurtis　距離
        #dMatrix[i, j] = distance.canberra(df.ix[i, 1:5], df.ix[j, 1:5]) # canberra　距離
        #dMatrix[i, j] = distance.chebyshev(df.ix[i, 1:5], df.ix[j, 1:5])  # chebyshev　距離
        #dMatrix[i, j] = distance.cityblock(df.ix[i, 1:5], df.ix[j, 1:5])  # cityblock　距離
        #dMatrix[i, j] = distance.correlation(df.ix[i, 1:5], df.ix[j, 1:5])  # correlation　距離
        #dMatrix[i, j] = distance.cosine(df.ix[i, 1:5], df.ix[j, 1:5])  # cosine　距離
        #dMatrix[i, j] = distance.euclidean(df.ix[i, 1:5], df.ix[j, 1:5])  # euclidean　距離
        #dMatrix[i, j] = distance.hamming(df.ix[i, 1:5], df.ix[j, 1:5])  # hamming　距離
        #dMatrix[i, j] = distance.jaccard(df.ix[i, 1:5], df.ix[j, 1:5])  # jaccard　距離
```


```python
#距離行列の中身を見る
pd.DataFrame(dMatrix).head()
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
      <th>...</th>
      <th>140</th>
      <th>141</th>
      <th>142</th>
      <th>143</th>
      <th>144</th>
      <th>145</th>
      <th>146</th>
      <th>147</th>
      <th>148</th>
      <th>149</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.000000</td>
      <td>0.035533</td>
      <td>0.040816</td>
      <td>0.051020</td>
      <td>0.009804</td>
      <td>0.055556</td>
      <td>0.035176</td>
      <td>0.014778</td>
      <td>0.068063</td>
      <td>0.040404</td>
      <td>...</td>
      <td>0.300000</td>
      <td>0.289855</td>
      <td>0.268482</td>
      <td>0.302817</td>
      <td>0.295775</td>
      <td>0.291971</td>
      <td>0.289575</td>
      <td>0.278810</td>
      <td>0.265455</td>
      <td>0.253846</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.035533</td>
      <td>0.000000</td>
      <td>0.026455</td>
      <td>0.026455</td>
      <td>0.035533</td>
      <td>0.090909</td>
      <td>0.041667</td>
      <td>0.030612</td>
      <td>0.032609</td>
      <td>0.015707</td>
      <td>...</td>
      <td>0.304029</td>
      <td>0.293680</td>
      <td>0.264000</td>
      <td>0.314079</td>
      <td>0.314079</td>
      <td>0.288390</td>
      <td>0.285714</td>
      <td>0.274809</td>
      <td>0.291045</td>
      <td>0.249012</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.040816</td>
      <td>0.026455</td>
      <td>0.000000</td>
      <td>0.021277</td>
      <td>0.040816</td>
      <td>0.096154</td>
      <td>0.026178</td>
      <td>0.035897</td>
      <td>0.038251</td>
      <td>0.031579</td>
      <td>...</td>
      <td>0.316176</td>
      <td>0.305970</td>
      <td>0.285141</td>
      <td>0.318841</td>
      <td>0.318841</td>
      <td>0.308271</td>
      <td>0.306773</td>
      <td>0.295019</td>
      <td>0.295880</td>
      <td>0.269841</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.051020</td>
      <td>0.026455</td>
      <td>0.021277</td>
      <td>0.000000</td>
      <td>0.051020</td>
      <td>0.096154</td>
      <td>0.026178</td>
      <td>0.035897</td>
      <td>0.027322</td>
      <td>0.021053</td>
      <td>...</td>
      <td>0.308824</td>
      <td>0.298507</td>
      <td>0.277108</td>
      <td>0.318841</td>
      <td>0.318841</td>
      <td>0.300752</td>
      <td>0.298805</td>
      <td>0.287356</td>
      <td>0.295880</td>
      <td>0.261905</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.009804</td>
      <td>0.035533</td>
      <td>0.040816</td>
      <td>0.051020</td>
      <td>0.000000</td>
      <td>0.055556</td>
      <td>0.035176</td>
      <td>0.014778</td>
      <td>0.068063</td>
      <td>0.040404</td>
      <td>...</td>
      <td>0.307143</td>
      <td>0.297101</td>
      <td>0.276265</td>
      <td>0.309859</td>
      <td>0.302817</td>
      <td>0.299270</td>
      <td>0.297297</td>
      <td>0.286245</td>
      <td>0.272727</td>
      <td>0.261538</td>
    </tr>
  </tbody>
</table>
<p>5 rows × 150 columns</p>
</div>



__距離行列を距離ベクトルに変換する__


```python
# 距離行列を距離ベクトルに変換する
dArray = distance.squareform(dMatrix)
```


```python
# 距離ベクトルを確認する
pd.DataFrame(dArray).T
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
      <th>...</th>
      <th>11165</th>
      <th>11166</th>
      <th>11167</th>
      <th>11168</th>
      <th>11169</th>
      <th>11170</th>
      <th>11171</th>
      <th>11172</th>
      <th>11173</th>
      <th>11174</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.035533</td>
      <td>0.040816</td>
      <td>0.05102</td>
      <td>0.009804</td>
      <td>0.055556</td>
      <td>0.035176</td>
      <td>0.014778</td>
      <td>0.068063</td>
      <td>0.040404</td>
      <td>0.028571</td>
      <td>...</td>
      <td>0.045593</td>
      <td>0.014749</td>
      <td>0.031884</td>
      <td>0.042424</td>
      <td>0.030864</td>
      <td>0.054545</td>
      <td>0.034921</td>
      <td>0.035294</td>
      <td>0.027692</td>
      <td>0.045317</td>
    </tr>
  </tbody>
</table>
<p>1 rows × 11175 columns</p>
</div>



__階層的クラスタリング__


```python
#クラスタリングする。リンケージ手法にもいろいろあるので、どれか選ぶ。
result = linkage(dArray, method = 'average')
#result = linkage(dArray, method = 'complete')
#result = linkage(dArray, method = 'single')
#result = linkage(dArray, method = 'weighted')
```

__デンドログラムの描画__


```python
# デンドログラムの描画
dend = dendrogram(result)
```


![output_29_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/f3cb8949-12d4-fef6-ce70-836c00a0abb9.png)


----------
## 課題
以下の課題を解いてください。
* __課題1__:  上記のクラスタリングで、「距離の定義」(metric)と「リンケージ手法」(method)を変えるとデンドログラムがどのように変わるか確認してください。
* __課題2__: 同様に、クラスタ数を２〜６に変えるとマッピング結果がどのように変わるか確認してください。