mnistをvariational auto encoderで教師なし学習し、潜在変数をクラスタリングして評価する

# mnistをvariational auto encoderで教師なし学習し、最終段をクラスタリングして評価する

```python
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
# 必要なライブラリのインポート
from keras.datasets import mnist
import numpy as np
import pandas as pd
import sklearn
# Jupyter notebookを利用している際に、notebook内にplot結果を表示するようにする
import matplotlib.pyplot as plt
%matplotlib inline

from keras.layers import Lambda, Input, Dense
from keras.models import Model
from keras.losses import mse
from keras import backend as K
import gc
```

    Using TensorFlow backend.



```python
feature_dims = range(2, 12)
```


```python
#Kerasの関数でデータの読み込み。データをシャッフルして学習データと訓練データに分割
(x_train, y_train), (x_test, y_test) = mnist.load_data()

# 2次元データを数値に変換
x_train = x_train.reshape(60000, 784)
x_test = x_test.reshape(10000, 784)
# 型変換
x_train = x_train.astype('float32')
x_test = x_test.astype('float32')
# 255で割ったものを新たに変数とする
x_train /= 255
x_test /= 255

# one-hot encodingを施すためのメソッド
from keras.utils.np_utils import to_categorical
# クラス数は10
num_classes = 10
y_train = y_train.astype('int32')
y_test = y_test.astype('int32')
labels = y_test
# one-hot encoding
y_train = to_categorical(y_train, num_classes)
y_test =  to_categorical(y_test, num_classes)
```


```python
def fitting(feature_dim, x_train, y_train, x_test, y_test):
    original_dim = x_train.shape[1]

    input_shape = (original_dim, )
    latent_dim = feature_dim

    # Reparametrization Trick 
    def sampling(args):
        z_mean, z_logvar = args
        batch = K.shape(z_mean)[0]
        dim = K.int_shape(z_mean)[1]
        epsilon = K.random_normal(shape=(batch, dim), seed = 5) # ε
        return z_mean + K.exp(0.5 * z_logvar) * epsilon

    # VAEモデル構築
    inputs = Input(shape=input_shape)
    x1 = Dense(256, activation='relu')(inputs)  
    x2 = Dense(64, activation='relu')(x1) 
    z_mean = Dense(latent_dim)(x2)
    z_logvar = Dense(latent_dim)(x2)
    z = Lambda(sampling, output_shape=(latent_dim,))([z_mean, z_logvar])
    encoder = Model(inputs, [z_mean, z_logvar, z], name='encoder')
    encoder.summary()

    latent_inputs = Input(shape=(latent_dim,))
    x3 = Dense(64, activation='relu')(latent_inputs)  
    x4 = Dense(256, activation='relu')(x3)  
    outputs = Dense(original_dim, activation='sigmoid')(x4)
    decoder = Model(latent_inputs, outputs, name='decoder')
    decoder.summary()

    z_output = encoder(inputs)[2]
    outputs = [decoder(z_output),z_output]
    vae = Model(inputs, outputs, name='variational_autoencoder')

    # 損失関数
    # Kullback-Leibler Loss
    kl_loss = 1 + z_logvar - K.square(z_mean) - K.exp(z_logvar)
    kl_loss = K.sum(kl_loss, axis=-1)
    kl_loss *= -0.5
    # Reconstruction Loss
    reconstruction_loss = mse(inputs, outputs[0])
    reconstruction_loss *= original_dim

    vae_loss = K.mean(reconstruction_loss + kl_loss)
    vae.add_loss(vae_loss)
    vae.compile(optimizer='adam')
    vae.summary()
    history = vae.fit(x_train,
                    epochs=50,
                    batch_size=256,
                    validation_data=(x_test, None))
    result = vae.predict(x_test)

    K.clear_session() # ←これです
    gc.collect()
    from IPython.display import clear_output
    clear_output()
    return (history, vae, result)
```


```python
#model = fitting(10, x_train, y_train, x_test, y_test)
models = [None] * len(feature_dims)
histories = [None] * len(feature_dims)
dec_imgs = [None] * len(feature_dims)
results = [None] * len(feature_dims)
for i in range(len(feature_dims)):
    (histories[i], models[i], dec_imgs[i]) = fitting(feature_dims[i], x_train, y_train, x_test, y_test)
```


```python
for i in range(len(feature_dims)):
    print(feature_dims[i])
    result = dec_imgs[i]
    decoded_imgs = result[0]
    # テスト画像と変換画像の表示
    n = 10
    plt.figure(figsize=(10, 2))
    for j in range(n):
        # テスト画像を表示
        ax = plt.subplot(2, n, j+1)
        plt.imshow(x_test[i].reshape(28, 28))
        plt.gray()
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

        # 変換された画像を表示
        ax = plt.subplot(2, n, j+1+n)
        plt.imshow(decoded_imgs[i][0][j].reshape(28, 28))
        plt.gray()
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)
    plt.show()
    results[i] = result[1]
```


![output_5_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/4411d80e-1685-67e6-d9ac-542a2f0b3a86.png)

    



![output_5_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/690e035f-8fa9-ad80-b0f7-51aa907ea14b.png)
    



![output_5_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/73cbccd4-59c0-3993-0a55-b4f5e8dd0eda.png)

    



![output_5_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/40877fa3-df57-89cb-12c2-4b4fd61a43f5.png)
    



![output_5_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/c0e47194-806c-9bcb-3fb4-7f005c31b1e4.png)
    



![output_5_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/0881c083-8fb3-7b00-324c-8a60d2f4a8f1.png)
    



![output_5_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/08bc7e5b-fe46-4f31-f07e-4b417f9014ba.png)

    


![output_5_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/da21f089-3509-805c-4c7d-8d0118374c18.png)
    



![output_5_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/b6d88f34-b3a9-d477-084c-59d1f7444c7d.png)

    



![output_5_9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/4a10248f-7712-7e00-79ec-0eca3e24a44a.png)
    



```python
for i in range(len(feature_dims)):
    results[i] = dec_imgs[i][1]
```


```python
#model.save('model/mnist-10')
#model = keras.models.load_model('model/mnist-10')
#for i in range(len(feature_dims)):
#    models[i].pop() # 最終段のsoftmax層を取り除いて、特徴量の層を最終段とする
#    models[i].summary()
```


```python
#result = model.predict(x_test)
#results = [None] * len(feature_dims)
#for i in range(len(feature_dims)):
#    keras.backend.clear_session()
#    results[i] = models[i].predict(x_test)
```


```python
def tsne(result):
    #t-SNEで次元削減
    from sklearn.manifold import TSNE
    tsne = TSNE(n_components=2, random_state = 0, perplexity = 30, n_iter = 1000)
    return tsne.fit_transform(result)
```


```python
#tsne = tsne(result)
tsnes = [None] * len(feature_dims)
for i in range(len(feature_dims)):
    tsnes[i] = tsne(results[i])
```


```python
#df = pd.DataFrame(tsne, columns = ['x', 'y'])
#df['label'] = labels
```


```python
def km(n_clusters, result):
    # k-meansでクラスタリングする
    from sklearn.cluster import KMeans
    return KMeans(n_clusters).fit_predict(result)
```


```python
#km = km(10, result)
#df['km'] = km
kms = [None] * len(feature_dims)
for i in range(len(feature_dims)):
    kms[i] = km(10, results[i])
```


```python
def DBSCAN(n_clusters, result):
    from sklearn.cluster import DBSCAN
    db = DBSCAN(eps=0.2, min_samples=n_clusters).fit(result)
    return db.labels_
```


```python
#dbscan = DBSCAN(20, result)
#df['DBSCAN'] = dbscan
```


```python
def hierarchy(result):
    from scipy.cluster.hierarchy import linkage, dendrogram
    result1 = linkage(result, 
                  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')
    return result1
```


```python
#hierarchy = hierarchy(result)
#display(hierarchy)
```


```python
def label_to_colors(label):
    color_dict = dict([(color[0], color[1]['color']) for color in zip(np.unique(label), plt.rcParams['axes.prop_cycle'])])
    colors = np.empty(label.shape, np.object)
    for k, v in color_dict.items():
        colors[label==k] = v
    return colors

#def cluster_visualization(x, y, label, cluster, method, n_clusters):
def cluster_visualization(x, y, label, cluster):
    plt.figure(figsize = (30, 15))
    plt.subplot(1,2,1)
    plt.scatter(x, y, c=label_to_colors(label))
#    for i in range(10):
#        tmp_df = df[df['label'] == i]
#        plt.scatter(tmp_df['x'], tmp_df['y'], label=i)
#    plt.legend(loc='upper left', bbox_to_anchor=(1,1))
    plt.subplot(1,2,2)
    plt.scatter(x, y, c=label_to_colors(cluster))
#    for i in range(n_clusters):
#        tmp_df = df[df[method] == i]
#        plt.scatter(tmp_df['x'], tmp_df['y'], label=i)
#    plt.legend(loc='upper left', bbox_to_anchor=(1,1))
```


```python
for i in range(len(feature_dims)):
    cluster_visualization(tsnes[i][:,0], tsnes[i][:,1], labels, kms[i])
```


![output_19_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/493adb16-9529-4a17-8070-0b224400b318.png)




![output_19_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/6a2e3d68-73ba-f288-e12f-d51dee999fdb.png)
    



![output_19_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/3e5c9091-9c2e-625d-f10b-0fc1ebfb62cf.png)
    



![output_19_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/23a3455e-c532-ccc4-8fc5-6443f6081f97.png)
    



![output_19_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/b87e10e2-93bc-13d5-add5-425024a40012.png)
    



![output_19_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/52fd8b1e-4d43-4f6e-cd76-21295ddb238b.png)
    



![output_19_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/9c53b83d-7d67-dacf-647f-ce632138bdd1.png)
    



![output_19_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/fccead80-1827-4a57-a127-a73099c07c77.png)
    



![output_19_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/32811362-5aa4-7886-b92d-ea2d4794df8a.png)
    



![output_19_9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/05e2c30f-e1a3-2bca-d398-278186597ab5.png)
    



```python
# https://qiita.com/mamika311/items/75c24f6892f85593f7e7
from sklearn.metrics.cluster import adjusted_rand_score
for i in range(len(feature_dims)):
    print("dim:" + str(feature_dims[i]) + " ARI: " + str(adjusted_rand_score(labels, kms[i])))
```

    dim:2 ARI: 0.36620498031529986
    dim:3 ARI: 0.41914836520424376
    dim:4 ARI: 0.49394921137719777
    dim:5 ARI: 0.5245649990462847
    dim:6 ARI: 0.47705674510916446
    dim:7 ARI: 0.41013993209378174
    dim:8 ARI: 0.3698302406743967
    dim:9 ARI: 0.32840225806718926
    dim:10 ARI: 0.4466812382927318
    dim:11 ARI: 0.4090677997413063



```python
# https://scikit-learn.org/stable/modules/generated/sklearn.metrics.normalized_mutual_info_score.html
# https://qiita.com/kotap15/items/38289edfe822005e1e44
from sklearn.metrics import normalized_mutual_info_score
#display(normalized_mutual_info_score(labels, df['km']))
for i in range(len(feature_dims)):
    print("dim:" + str(feature_dims[i]) + " NMI: " + str(normalized_mutual_info_score(labels, kms[i])))
```

    dim:2 NMI: 0.5199419992579754
    dim:3 NMI: 0.56100353575167
    dim:4 NMI: 0.605060303081276
    dim:5 NMI: 0.6020900415664949
    dim:6 NMI: 0.5631744057166579
    dim:7 NMI: 0.5014462787979749
    dim:8 NMI: 0.46110014862882315
    dim:9 NMI: 0.42836636346088663
    dim:10 NMI: 0.5187118150024308
    dim:11 NMI: 0.48519256224162205



```python
def shilhouette(clusters, x_test):
    from sklearn.metrics import silhouette_samples
    from matplotlib import cm
    plt.figure(figsize = (10, 10))
    cluster_labels=np.unique(clusters)
    n_clusters=cluster_labels.shape[0]
    silhouette_vals=silhouette_samples(x_test,clusters,metric='euclidean')
    y_ax_lower,y_ax_upper=0,0
    yticks=[]
    for i,c in enumerate(cluster_labels):
        c_silhouette_vals=silhouette_vals[clusters==c]
        print(len(c_silhouette_vals))
        c_silhouette_vals.sort()
        y_ax_upper +=len(c_silhouette_vals)
        color=cm.jet(float(i)/n_clusters)
        plt.barh(range(y_ax_lower,y_ax_upper),
                c_silhouette_vals,
                height=1.0,
                edgecolor='none',
                color=color
                )
        yticks.append((y_ax_lower+y_ax_upper)/2.)
        y_ax_lower += len(c_silhouette_vals)

    #シルエット係数が1であれば　良くクラスタリングできてる 
    #またシルエットの幅がクラスタ数で平均して等しいとき、全体のデータを等分割できていることを示す
    #この分割幅=シルエットバーの幅が等しくなり、かつ、シルエット係数が1に近づくようにkを最適化することが設定手法として考えられる.

    #平均の位置に線を引く
    silhouette_avg=np.mean(silhouette_vals)
    plt.axvline(silhouette_avg,color="red",linestyle="--")
    plt.ylabel("Cluster")
    plt.xlabel("Silhouette coefficient")
```


```python
for i in range(len(feature_dims)):
    shilhouette(kms[i], x_test)
```

    1391
    704
    1316
    843
    582
    1102
    500
    1330
    848
    1384
    1196
    1035
    1287
    901
    976
    751
    718
    859
    695
    1582
    852
    657
    1130
    1615
    977
    994
    904
    1186
    914
    771
    883
    1316
    1130
    964
    982
    928
    737
    829
    1140
    1091
    997
    898
    1000
    1122
    911
    1078
    803
    1152
    998
    1041
    1149
    1211
    976
    1106
    765
    1100
    852
    1010
    886
    945
    1070
    696
    1395
    825
    934
    1008
    926
    969
    981
    1196
    1117
    1060
    975
    1118
    807
    970
    1007
    922
    881
    1143
    921
    878
    948
    971
    1033
    1098
    1109
    1077
    1013
    952
    1008
    888
    1248
    1059
    1166
    848
    925
    961
    889
    1008



    
![output_23_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/06f32dfc-dd82-2b29-d5a0-ed2ad8adf28a.png)
    



![output_23_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/103977b8-253d-e9b0-0744-44e94737c0ff.png)
    



![output_23_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/eaf83723-0670-3952-2a2a-3237cca48745.png)
    



![output_23_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/9637761b-311d-5f5b-dd18-5a1eacbbfc72.png)
    



![output_23_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/feed1c87-c491-92e8-6c8f-cf8846193e88.png)
    



![output_23_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/fbd0ea92-bfaf-8d79-0810-6558f191280e.png)
    



![output_23_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/c9666481-f15e-00ce-59e3-eff7edbea338.png)
    



![output_23_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/cf26a845-9ba9-fc04-f295-45b065e21a7f.png)
    



![output_23_9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/be56e031-64c3-cc35-e855-9ce41eac0b6f.png)

    



![output_23_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/490869/d17da16f-60be-8469-0554-dc0f82ff02d6.png)
    



```python

```