1:ライブラリのインポート

今回の演習で用いる基本的なライブラリを以下のようインポートし、使用可能な状態にする。

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sklearn
import time

cmap='tab10' #グラフをプロットしたときのカラーコードの指定

numpyは数値計算を効率的に行うためのライブラリであり、

impport numpy as np

は、numpyというライブラリをインポートし(使用可能とし)、それを以降「np」という略号で参照するという意味になる。

同様に、matplotlibはグラフのプロットに用いるライブラリ、pandasもデータ解析時に便利なコマンドを多く含むライブラリである。また、sklearnは、初心者向けの機械学習のオープンソースライブラリである。

上でインポートしたのは共通で用いる基本的なライブラリであり、特定の解析で用いるライブラリは適宜そこでインポートする。

2:データのダウンロード

今回の演習では、最初に手書き文字の画像データであるMNISTを用い、データの読み込みとそこからの解析を行う。MINSTは28×28ピクセルの画像ファイルであり、機械学習の練習台として広く用いられる。googleで検索すると沢山ヒットするが、例えば以下のサイトなどを参照のこと。

MNISTって何?

まず、以下のリンクからmnist_train.csvとmnist_test.csvのファイルをダウンロードして保存する。

mnist_train.csv

mnist_test.csv

ダウンロードしたファイルは、notebookを開いたフォルダに移動させておくと良い。

Excelなどで開いてみるとわかるが、このcsvファイルは最初の1列に0~9の数字があり、これが正解ラベルとなっている。その後の2行目から785行目は、28×28=784個のピクセル強度の値であり、0~255の値をとっている。ファイルは60001行からなり、一番上の行はヘッダーとして列の説明のために用いるものであり、それを取り除いた60000個のデータを含んでいる。

以降、mnist_trainのデータを用いて主成分分析などの解析を行う。mnist_testは予測精度の評価などに用いる。

3:データの読み込み

では、MNIST_train.csvのデータを読み込んでみる。ここではpandasライブラリにあるread_csvという関数を用いる。

In [2]:
xx = np.array(pd.read_csv('./mnist_train.csv'))

カッコの中はファイル名である。その前にある「./」はこのNotebookがあるフォルダという意味を持つ(相対パスと呼ぶ)。異なるフォルダにおいたファイルを読み込む場合には、以下のように絶対パスを用いる。

np.array関数でくくっているのは、読み込んだデータをnumpyのarray形式に変換するためであるが、あまり気にしなくて良い。

In [3]:
xx = np.array(pd.read_csv('/Users/chkar/Desktop/VAE/mnist_train.csv'))

もちろん、上のフォルダの位置は環境に依存するので、各自適切なところを指定する。

上のread_csv関数は、デフォルトでヘッダーの部分として一番上の行を取り除く。なので上の変数xxは、データ数の60000行×785列のデータを含む。それを確認するには、

In [4]:
xx.shape
Out[4]:
(60000, 785)

と変数の後にピリオドを付けて「shape」と打てば良い(これは、np.shape(xx)と同じ意味になる)。

次にこの変数xxから、ラベル情報(1列目)を抜き出す操作を行う。このとき、画像の縦横のピクセル数の28は頻出する値なので、わかりやすい変数名を付けて定義をしておくことにする

In [5]:
pixel_size=28

x_train = xx[:,1:].astype('float32') / 255
y_train = xx[:,0] 

これでx_trainは、csvファイルの行列から最初の1列を取り除いた60000行x784列のデータになる。この配列の操作は、2018年12月の理論情報交換会の資料

Python入門

を参考にするとよい。上のxx[:,1:]にある「1:」は、1列目以降全てという意味になる(Pythonでは最初の列は0列目になることに注意)。

上のコマンド群により、2次元配列x_trainに画像データが、1次元配列y_trainに正解ラベルが入ることになる(x_train.shapeなどにより配列のサイズを確認すると良い)。「.astype('float32') / 255」では、元々のデータが0~255の整数なので、それを0~1の値に変換をしている。

データの読み込みについて、上のようにコマンドを繋げて処理をしても良いし、あるいは下のように読み込むための関数を自分で定義して、それを用いても良い。

In [6]:
def load_csv(csv):
    xx=np.array(pd.read_csv(csv))
    x_data=xx[:,1:].astype('float32') / 255
    y_data=xx[:,0]
    
    return x_data, y_data

この関数を定義することにより、以下のようにデータを読み込むことが出来る。

In [7]:
x_train, y_train=load_csv('./mnist_train.csv')
x_test, y_test=load_csv('./mnist_test.csv')

こうしておくと、このload_csv関数を後で別のデータの読み込みに使い回すことが出来る。

このMNISTのデータはtrainが60000個、testが10000個のデータを含む(.shapeを付けて確認してみると良い)。このデータ数だと、後から行う解析のいくつかは(今回のハンズオンセミナーでは)実行時間が長くなり過ぎる。なので、上のファイルからデータ数を1/10にしたデータセットを以下のように作成する。このサイズのデータセットでも、自分のPCでは遅いなと感じたら、1/20などより少ない数のデータからなるデータセットを作成する。

In [8]:
n=int(x_train.shape[0]/10) #データの個数(x_train.shape[0])を10で割って整数にしている
x_train2=x_train[:n,]
y_train2=y_train[:n]

n=int(x_test.shape[0]/10) 
x_test2=x_test[:n,]
y_test2=y_test[:n]

今回の演習では、MNISTの場合には処理時間の関係から主にデータ数を6000にしたx_train2, y_train2を用る。

また、MINISTとは違う画像データセットとして、以下の2つのデータセットを用意した。

Fashion MNIST

植物葉の形態データ

上のリンクからダウンロードできる。それぞれの解析について、まずMNISTで試した後に、上の2つのデータセットでも試してみるのが良い。ただし、MNIST, Fashion MNISTの画像は28x28 pixelであるのに対して、植物葉の形態データは64x64 pixelであり、適宜パラメータを変える必要がある(上で定義したpixel_sizeの値を必要に応じて変更する)。

4:画像のチェック

読み込んだ画像データがどのようなものかチェックをしてみる。以下の関数で、3次元配列であるx_mnist_trainの0番目のデータを画像として出力することが出来る。下で使う関数plt.inshowで画像を表示できるように、x_trainの784要素のベクトルを、28x28の形に以下のように変換する。

In [9]:
x_img=x_train.reshape([x_train.shape[0],pixel_size,pixel_size])

これでx_train2は60000x28x28の大きさの3次元配列となる。その最初のデータを以下のように表示してみる。

In [10]:
plt.imshow(x_img[0,:], cmap='gray')
Out[10]:
<matplotlib.image.AxesImage at 0x1c4237a91c8>

複数の画像を一度に見たいときには、以下のようにすると良い。最初のplt.figure関数で新規のウィンドウを描画し、その後にforループを用いてadd_subplot関数により図を加えている。

In [11]:
fig = plt.figure(figsize=(12,10))
for i in range(0,25):
    axes = fig.add_subplot(5, 5, i%25+1)
    axes.imshow(x_img[i,:,:],cmap='gray')

上のセルを実行してみると、手書きの数字ではあるが、なかなか汚いことが見て取れる。

5:主成分分析

では、上の画像データセットを主成分分析(Principle Component Analysis; PCA)により解析してみる。まず、用いるライブラリを導入する。

In [12]:
from sklearn.decomposition import PCA

以下の関数群によりPCAを実行する

In [13]:
pca=PCA(n_components=8)
transformed=pca.fit_transform(x_train2)

「pca=PCA(n_components=8)」は、PCAクラスのインスタンスを生成しているのだが、あまり深くは追求はしない。クラスは、ざっくりと言うとデータや処理の定義をまとめた雛型であり、インスタンスはそのクラスを実体化したもので、クラスで定義された変数や関数を操作することが可能となる。また、カッコ内のn_componentは何番目の主成分まで見るかを決める変数である。

fit_transform関数がPCAを実施する関数であり、その結果がdata_size×n_componentの大きさを持つ配列として出力される。

結果をプロットするには、以下のようにすると良い。

In [14]:
plt.scatter(transformed[:,0],transformed[:,1],s=3,c=y_train2,cmap=cmap)
plt.colorbar()
plt.xlabel("PC1")
plt.ylabel("PC2")
Out[14]:
Text(0, 0.5, 'PC2')

上のセルにおいて、関数plt.scatterの引数であるtransformed[:0,]とtransformed[:1,]はそれぞれ第一主成分と第二主成分である。それがグラフのx軸とy軸のデータになる。以降の引数の「s」は点のサイズ、「c=y_train2」はデータ点のカテゴリがベクトルy_train2で与えられることを意味する。またcmapでカラーマップ 「tab10」を指定している。

(カラーバーの色と値が多少ずれますが、無視してください。)

どうせなら、3次元の主成分空間でどのように点が分布をしているかを見たいかもしれない。その場合、プロットをぐりぐり回せるのが良い。その方法については、以下のリンクを参照。

3D plot

6:t-SNEによる可視化

t-SNE は、高次元データの次元を圧縮するアルゴリズムであり、Single-cell RNA-Seqなど高次元のデータの可視化によく用いられる。手法の詳細は講義資料を参考のこと。ここでは、t-SNEを用いた画像データの次元圧縮を行ってみる。

まずは必要なライブラリをインポートする(sklearnのライブラリを用いる)。

In [15]:
from sklearn.manifold import TSNE

次に下の関数でクラスTSNEをインスタンス化する(意味は深く考える必要なない)。n_componentsで何次元まで圧縮するかを指定している。

In [16]:
tsne = TSNE(n_components=2)

t-SNEの計算にはすこし時間がかかる。1/10にしたデータセット(x_train2, y_train2)を使うのが無難である。あるいはもっと少ないデータ数まで減らしても良い。

下のセルでは、2行目がt-SNEの計算をする関数であるが、その前後でtime関数を用いてその時刻を記録し、最後に計算に要した時間を表示している(単位は秒)。

In [17]:
t1 = time.time() 
transformed = tsne.fit_transform(x_train2)
t2 = time.time()
print(f"経過時間:{t2-t1}")
経過時間:87.11753273010254

古澤が所有するThinkpad X1 carbon(2018年モデル/core i5-8250U)で、1/10にしたMINSTデータの処理に90秒程度を要している。

下のセルの関数で、2次元平面にt-SNEによって圧縮した結果をプロットできる。

In [18]:
plt.scatter(transformed[:,0],transformed[:,1],s=3,c=y_train2,cmap=cmap)
plt.colorbar()
Out[18]:
<matplotlib.colorbar.Colorbar at 0x1c441bf76c8>

発展課題 6-1:
「n_components=3」とすることにより、3次元空間にt-SNEによりデータを次元圧縮することが出来る。上の3d_plotリンクを参考にして、3次元空間でどのようにデータが分布をするか調べてみよ。

発展課題 6-2:
t-SNEにはperplexityという重要なパラメータがある。これは、ざっくりと言うとどれくらい近傍の点を計算に取り入れるか、ということを決めているパラメータである。この値は、

tsne = TSNE(n_components=2, perplexity=30)

といった形で指定をすることが出来る(デフォルトは30)。この値を、1,3, 30, 100などと1-100の範囲で変化をさせてみて、結果がどのように変わるかを観察してみよ。

発展課題 6-3
t-SNEの弱点の一つは、計算時間が比較的長い点にある。計算時間を短くする一つの方法として、まずデータをPCAにかけて、その上位の主成分だけを用いてt-SNEを用いる場合がある。つまりMNISTの場合は元々は28x28=784次元あるが、その全てが情報を持っているわけではない(例えば、隅っこのデータはいつでも黒とか)。なので、PCAを用いて上位の10次元も取れば、十分にt-SNEで分離できる場合がある。それを試してみて、結果のグラフがどのように変化するか、計算時間がどれくらい削減できるかを調べてみよ。

7:UMAPによる可視化

UMAPは2018年に発表された比較的新しい手法であり、t-SNEと同じような次元圧縮を数倍の速さで実行できる。おそらく、今後t-SNEに変わり用いられていく手法になると思われる。まずはライブラリんをインポートする。

In [19]:
import umap
import warnings
warnings.filterwarnings('ignore') #warningが出るので無視するようにしてしまっている

from scipy.sparse.csgraph import connected_components

以下のようにUMAPを実行する。

In [20]:
t1 = time.time() 
result = umap.UMAP(n_neighbors=5, n_components=2).fit(x_train2)
t2 = time.time()
print(f"経過時間:{t2-t1}")
経過時間:18.910237789154053

上の操作により、resultに結果が格納される。そこから2次元に圧縮されたデータを取り出すには、result.embedding_という配列を以下のようにプロットする。

In [21]:
plt.figure(figsize=(6,5))
plt.scatter(result.embedding_[:,0],result.embedding_[:,1],s=3,c=y_train2, cmap=cmap)
plt.colorbar()
Out[21]:
<matplotlib.colorbar.Colorbar at 0x1c44f9e8a88>

ある程度は綺麗に分離できていることがわかるはずである。

発展問題 7-1:
t-SNEの場合と同じように、「n_components=3」とすることにより、3次元空間にデータを圧縮することが出来る。上の3d_plotリンクを参考にして、3次元空間でどのようにデータが分布をするか調べてみよ。

発展問題 7-2:
UMAPの利点の一つは、次元圧縮した後の空間に、新たなデータ点を加えることが出来る点にある(t-SNEは原理的にそれが出来ない)。トレーニングデータであるx_train2から作られた空間に、テストデータであるx_test2をプロットするには、以下のようにする。

test_embedding=result.transform(x_test2)

そして、x_test2のデータがどのように分布するかを見たければ、test_embedding[:,0]とtest_embedding[:,1]をplt.scatter関数でプロットすれば良い。同じような分離が出来ているか、チェックしてみよ。

8:変分オートエンコーダ(VAE)による解析

まず、必要なライブラリを以下のようにインポートする。

In [22]:
from keras.layers import Lambda, Input, Dense
from keras.models import Model
from keras.losses import mse, binary_crossentropy
from keras import backend as K
Using TensorFlow backend.

VAEの解析では、関数群が複雑になるので、以下のように自分で関数を新たに定義して解析を進める。これは、Tensorflow/KerasライブラリにあるVAEのサンプルコードを適宜改変したものである。

In [23]:
def vae_mlp(x_train, x_test, latent_dim = 2, epochs = 30, pixel_size=pixel_size):

    # network parameters
    input_shape = (pixel_size*pixel_size, )
    intermediate_dim = 512
    batch_size = 128

    # build encoder model
    inputs = Input(shape=input_shape, name='encoder_input')
    x = Dense(intermediate_dim, activation='relu')(inputs)
    z_mean = Dense(latent_dim, name='z_mean')(x)
    z_log_var = Dense(latent_dim, name='z_log_var')(x)

    # use reparameterization trick to push the sampling out as input
    z = Lambda(sampling, output_shape=(latent_dim,), name='z')([z_mean, z_log_var])

    # instantiate encoder model
    encoder = Model(inputs, [z_mean, z_log_var, z], name='encoder')
    encoder.summary()

    # build decoder model
    latent_inputs = Input(shape=(latent_dim,), name='z_sampling')
    x = Dense(intermediate_dim, activation='relu')(latent_inputs)
    outputs = Dense(pixel_size*pixel_size, activation='sigmoid')(x)

    # instantiate decoder model
    decoder = Model(latent_inputs, outputs, name='decoder')
    decoder.summary()

    # instantiate VAE model
    outputs = decoder(encoder(inputs)[2])
    vae = Model(inputs, outputs)

    # VAE loss = mse_loss or xent_loss + kl_loss
    reconstruction_loss = binary_crossentropy(inputs,outputs)*pixel_size*pixel_size

    kl_loss = 1 + z_log_var - K.square(z_mean) - K.exp(z_log_var)
    kl_loss = -0.5*K.sum(kl_loss, axis=-1)
    
    vae_loss = K.mean(reconstruction_loss + kl_loss)
    vae.add_loss(vae_loss)
    vae.compile(optimizer='adam')
    vae.summary()

    result=vae.fit(x_train,
             epochs=epochs,
             batch_size=batch_size,
             validation_data=(x_test, None))
    
    return encoder, decoder, result

def sampling(args):
    z_mean, z_log_var = args
    batch = K.shape(z_mean)[0]
    dim = K.int_shape(z_mean)[1]
    epsilon = K.random_normal(shape=(batch, dim))
    return z_mean + K.exp(0.5 * z_log_var) * epsilon

この関数vae_mlpを用いて、x_train2のデータを以下のように学習させる。

In [24]:
encoder, decoder, result= vae_mlp(x_train2, x_test2, latent_dim = 2, epochs = 30, pixel_size=pixel_size)
Model: "encoder"
__________________________________________________________________________________________________
Layer (type)                    Output Shape         Param #     Connected to                     
==================================================================================================
encoder_input (InputLayer)      (None, 784)          0                                            
__________________________________________________________________________________________________
dense_1 (Dense)                 (None, 512)          401920      encoder_input[0][0]              
__________________________________________________________________________________________________
z_mean (Dense)                  (None, 2)            1026        dense_1[0][0]                    
__________________________________________________________________________________________________
z_log_var (Dense)               (None, 2)            1026        dense_1[0][0]                    
__________________________________________________________________________________________________
z (Lambda)                      (None, 2)            0           z_mean[0][0]                     
                                                                 z_log_var[0][0]                  
==================================================================================================
Total params: 403,972
Trainable params: 403,972
Non-trainable params: 0
__________________________________________________________________________________________________
Model: "decoder"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
z_sampling (InputLayer)      (None, 2)                 0         
_________________________________________________________________
dense_2 (Dense)              (None, 512)               1536      
_________________________________________________________________
dense_3 (Dense)              (None, 784)               402192    
=================================================================
Total params: 403,728
Trainable params: 403,728
Non-trainable params: 0
_________________________________________________________________
Model: "model_1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
encoder_input (InputLayer)   (None, 784)               0         
_________________________________________________________________
encoder (Model)              [(None, 2), (None, 2), (N 403972    
_________________________________________________________________
decoder (Model)              (None, 784)               403728    
=================================================================
Total params: 807,700
Trainable params: 807,700
Non-trainable params: 0
_________________________________________________________________
Train on 6000 samples, validate on 1000 samples
Epoch 1/30
6000/6000 [==============================] - 1s 191us/step - loss: 309.0391 - val_loss: 219.1690
Epoch 2/30
6000/6000 [==============================] - 1s 129us/step - loss: 212.8271 - val_loss: 200.2263
Epoch 3/30
6000/6000 [==============================] - 1s 134us/step - loss: 198.4307 - val_loss: 191.2396
Epoch 4/30
6000/6000 [==============================] - 1s 140us/step - loss: 190.6261 - val_loss: 186.0994
Epoch 5/30
6000/6000 [==============================] - 1s 141us/step - loss: 185.8243 - val_loss: 182.4785
Epoch 6/30
6000/6000 [==============================] - 1s 136us/step - loss: 182.3750 - val_loss: 179.4839
Epoch 7/30
6000/6000 [==============================] - 1s 159us/step - loss: 179.5138 - val_loss: 177.6498
Epoch 8/30
6000/6000 [==============================] - 1s 134us/step - loss: 176.5036 - val_loss: 174.6445
Epoch 9/30
6000/6000 [==============================] - 1s 134us/step - loss: 173.8697 - val_loss: 172.8189
Epoch 10/30
6000/6000 [==============================] - 1s 134us/step - loss: 172.0518 - val_loss: 172.4589
Epoch 11/30
6000/6000 [==============================] - 1s 128us/step - loss: 170.7633 - val_loss: 170.9805
Epoch 12/30
6000/6000 [==============================] - 1s 127us/step - loss: 169.7795 - val_loss: 170.5333
Epoch 13/30
6000/6000 [==============================] - 1s 126us/step - loss: 169.0404 - val_loss: 169.5901
Epoch 14/30
6000/6000 [==============================] - 1s 128us/step - loss: 168.3394 - val_loss: 169.3128
Epoch 15/30
6000/6000 [==============================] - 1s 125us/step - loss: 167.7386 - val_loss: 169.0956
Epoch 16/30
6000/6000 [==============================] - 1s 131us/step - loss: 167.3650 - val_loss: 168.7151
Epoch 17/30
6000/6000 [==============================] - 1s 131us/step - loss: 166.7812 - val_loss: 168.1580
Epoch 18/30
6000/6000 [==============================] - 1s 148us/step - loss: 166.3274 - val_loss: 168.7324
Epoch 19/30
6000/6000 [==============================] - 1s 141us/step - loss: 165.9151 - val_loss: 167.8224
Epoch 20/30
6000/6000 [==============================] - 1s 131us/step - loss: 165.4332 - val_loss: 166.9874
Epoch 21/30
6000/6000 [==============================] - 1s 126us/step - loss: 165.0775 - val_loss: 166.7476
Epoch 22/30
6000/6000 [==============================] - 1s 127us/step - loss: 164.6218 - val_loss: 167.1901
Epoch 23/30
6000/6000 [==============================] - 1s 126us/step - loss: 164.2430 - val_loss: 166.2464
Epoch 24/30
6000/6000 [==============================] - 1s 127us/step - loss: 163.7277 - val_loss: 165.9477
Epoch 25/30
6000/6000 [==============================] - 1s 126us/step - loss: 163.4324 - val_loss: 165.6998
Epoch 26/30
6000/6000 [==============================] - 1s 125us/step - loss: 163.0690 - val_loss: 165.6006
Epoch 27/30
6000/6000 [==============================] - 1s 125us/step - loss: 162.6324 - val_loss: 165.4999
Epoch 28/30
6000/6000 [==============================] - 1s 128us/step - loss: 162.2901 - val_loss: 164.6601
Epoch 29/30
6000/6000 [==============================] - 1s 127us/step - loss: 161.9139 - val_loss: 164.3650
Epoch 30/30
6000/6000 [==============================] - 1s 128us/step - loss: 161.6712 - val_loss: 164.2422

このVAEによる学習がどのように進行したかを確認するためには損失関数(loss function)をプロットすると良い。損失関数は入力データと出力データの差を評価する関数であり、学習においてそれが最小となるようにネットワークの重みが決められる。上のモデルでは損失関数として交差エントロピーが使われてるが、これは入力と出力の確率分布の違いを定量したものである(と思えばよい)。

損失関数には、学習に用いたtraining dataに対するものと、そうでないtest dataに対するものがあり、それぞれをプロットするには以下のようにすると良い。

In [25]:
plt.plot(result.history['loss'],label='loss',color='r')
plt.plot(result.history['val_loss'],label='val_loss',color='b')
plt.legend()
plt.xlabel('epoch')
plt.ylabel('loss')
plt.show()

lossとあるのが、training dataに対する損失関数、val_lossとしたのがtest dataに対するものである。もしlossが下がりながらも、val_lossが上昇している場合には、いわゆる過学習が生じている状況となる。

では、この学習したencoderを用いて、x_train2のデータがどのような空間にマップされているかを描いてみる。以下の命令で、配列z_meanにそれぞれのデータが2次元の潜在空間(Latent space)のどこにマップされるかが格納される(アンダーバーは、その値は無視して捨てることを意味する)。

In [26]:
z_mean, _, _ = encoder.predict(x_train2)

plt.figure(figsize=(6, 5))
plt.scatter(z_mean[:, 0], z_mean[:, 1], s=3, c=y_train2, cmap=cmap)
plt.colorbar()
Out[26]:
<matplotlib.colorbar.Colorbar at 0x1c44c2845c8>

次に、学習したdecoderを用いて、2次元の潜在空間でそれぞれの図形がどのように表現されているかを見る。そのために学習したdecoderを入力した以下の関数を用意する。この関数では、潜在空間の座標を-4から4まで変化させたときに、どのような形が出力されるかを描くものである。

In [27]:
def plot_latent(decoder):
    
    n=30
    figure = np.zeros((pixel_size * n, pixel_size * n))
 
    grid_x = np.linspace(-4, 4, n)
    grid_y = np.linspace(-4, 4, n)[::-1]    
  

    for i, yi in enumerate(grid_y):
        for j, xi in enumerate(grid_x):
            z_sample = np.array([[xi, yi]])
            x_decoded = decoder.predict(z_sample)
            digit = x_decoded[0].reshape(pixel_size, pixel_size)
            figure[i * pixel_size: (i + 1) * pixel_size,
                   j * pixel_size: (j + 1) * pixel_size] = digit

    plt.figure(figsize=(10, 10))
    start_range = pixel_size // 2
    end_range = (n - 1) * pixel_size + start_range + 1
    pixel_range = np.arange(start_range, end_range, pixel_size)
    sample_range_x = np.round(grid_x, 1)
    sample_range_y = np.round(grid_y, 1)
    plt.xticks(pixel_range, sample_range_x)
    plt.yticks(pixel_range, sample_range_y)
    plt.xlabel("z[0]")
    plt.ylabel("z[1]")
    plt.imshow(figure, cmap='Greys_r')
    plt.show()
    
plot_latent(decoder)    

発展問題 8-1:
上のVAEの解析は、データ数を1/10に削減したx_train2, y_train2のデータセットでは学習のためのデータ量が足りていない。削減前のx_train, y_trainで試すと、よりクリアにlatent spaceと数字画像との対応を見出すことができるようである。試してみよ(10分以上かかるかもしれない)。

発展問題 8-2:
上のvae_mlp関数では、encoder, decoderともに非常に単純なニューラルネットワークを用いてる(784次元のベクトルを512次元のベクトルに変換し(ここで784x512=40148個の重み係数が必要となる)、そこからz_meanとz_log_varを出力している)。より気の利いた、画像の解析を得意とするニューラルネットとして、下のsection 11で説明するConvolutional Neural Network(CNN)。がある。このCNNをVAEのencoder, decoderとして用いたプログラムを、

VAE with CNN

に置くので、試してみよ。その差が明確に見えるのは、おそらく1/10に削減しないデータを用いたときである。val_lossの落ち方や、再構成した画像の変化にも違いがあると思う。

9:Support Vector Machine(SVM)による分類

次に、SVMを用いてデータを分類する問題を考える。具体的には、テストデータを学習し、その学習に用いていないテストデータがどのカテゴリ(MNISTの場合は0~9の数字)に属するかを当てるという問題である。まず、ライブラリを以下のようにインポートする。

In [28]:
from sklearn import svm
from sklearn.metrics import confusion_matrix

SVMもけっこう時間がかかるので、データ数を1/10にしたx_train2を用いるのが良い。以下の命令を実行することにより、svm_clfというインスタンスに結果が格納される。

In [29]:
svm_clf = svm.SVC(gamma='auto')
svm_clf.fit(x_train2, y_train2)
Out[29]:
SVC(C=1.0, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma='auto', kernel='rbf',
    max_iter=-1, probability=False, random_state=None, shrinking=True,
    tol=0.001, verbose=False)

では、SVMでどれくらいの正解率が出たかをチェックしてみる。ここではaccuracyを用いる。これは、真陽性・偽陽性・真陰性・偽陰性をTP・FP・TN・FNとしたとき、

accuracy=(TP+TN)/(TP+FP+TN+FN)

で定義される量である。SVMの場合は、accuracyは以下のように出力できる。

In [30]:
svm_clf.score(x_test2, y_test2)
Out[30]:
0.885

次に、それぞれのカテゴリにおいて予測がどれくらい正しいかをVisualizeするために、以下の関数を用意してconfusion matrixを描くことにする。

In [31]:
def plot_confusion_matrix(cm,classes,cmap=plt.cm.Blues):
 
    plt.figure(figsize=(7,7))
    
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    
    plt.xticks(range(len(classes)),classes, rotation=45)
    plt.yticks(range(len(classes)),classes)
    plt.ylim(len(classes)-0.5, -0.5)
    
    thresh = cm.max() / 2.
    for i in range(cm.shape[0]): 
        for j in range(cm.shape[1]):
            plt.text(j, i, cm[i, j],
                     horizontalalignment="center",
                     color="white" if cm[i, j] > thresh else "black")
                 
    plt.colorbar(shrink=0.85)

    plt.ylabel('True label', fontsize=12)
    plt.xlabel('Predicted label' ,fontsize=12)

    plt.show()    

上の関数を用いてconfusion matrixを描くには以下のようにする。

In [32]:
prediction = svm_clf.predict(x_test2)
classifier = confusion_matrix(y_test2, prediction)

class_names = ['0','1','2','3','4','5','6','7','8','9']
plot_confusion_matrix(classifier,class_names)

発展問題 9-1:
SVMによる分類のaccuracyを改善する一つの方法として、適当な前処理によってデータを次元圧縮するという手法が考えられる。例えば、発展問題 7-2で行ったように、training dataをUMAPを用いて次元圧縮し、さらに同じ変換でtest dataも圧縮しておく。この圧縮後のデータを用いてSVMを行うと、accuracyが上昇する場合がある。MNISTの場合、UMAPのパラメータでn_neighbor=5, n_component=4くらいが良いかもしれない。試してみよ。

発展問題 9-2:
上の解析では、1/10にデータ数を削減したx_train2, x_test2を持ちて分類を行ったと思うが、削減前のx_trainを用いると、accuracyは上がる。つまり1/10ではSVMによる学習にはtraining dataが足りていない。興味があれば、削減前のx_trainのデータでSVMで分類してみよ。ただし、これは計算機によるが時間がかなり要することに注意すること(古澤のThinkpadで20分以上)。

10:Random Forestによる分類

次に、Random Forestによる分類を試してみる。必要なライブラリを以下のようにインポートする。

In [33]:
from sklearn.ensemble import RandomForestClassifier

SVMの時と同様に、以下の命令を実行することにより、randomforest_clfというインスタンスに結果が格納される。

In [34]:
randomforest_clf = RandomForestClassifier(n_estimators=10)
randomforest_clf.fit(x_train2, y_train2)
Out[34]:
RandomForestClassifier(bootstrap=True, ccp_alpha=0.0, class_weight=None,
                       criterion='gini', max_depth=None, max_features='auto',
                       max_leaf_nodes=None, max_samples=None,
                       min_impurity_decrease=0.0, min_impurity_split=None,
                       min_samples_leaf=1, min_samples_split=2,
                       min_weight_fraction_leaf=0.0, n_estimators=10,
                       n_jobs=None, oob_score=False, random_state=None,
                       verbose=0, warm_start=False)

実行に要する時間が、SVMよりもかなり短いはずである。accuracyを表示させるには、以下のようにする。

In [35]:
randomforest_clf.score(x_test2, y_test2)
Out[35]:
0.872
In [36]:
prediction = randomforest_clf.predict(x_test2)
classifier = confusion_matrix(y_test2, prediction)

class_names = ['0','1','2','3','4','5','6','7','8','9']
plot_confusion_matrix(classifier,class_names)

発展問題10-1
問題9-2と同様に、データ数の削減前のx_train, x_testを用いてRandom Forestでの分類を試してみよ。こちらは比較的短時間で終わると思われる。

11:Convolution Neural Network(CNN)による分類

では次に、CNNによる分類を試してみる。ここではKerasというTensoflow(Googleが開発している機械学習のライブラリ)を簡単に使えるようにするラッパーのライブラリを用いる。まずはライブラリのインポートを以下のように行う。

In [37]:
import keras
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras import backend as K

CNNだと、それぞれの入力は1次元のベクトルではなく、2次元配列(例えばMNISTの場合は28x28の画像ファイル)にしておく。これは、畳み込み(convolution)のときに2次元のフィルターを用いて処理をするなど、基本的に2次元画像の解析にライブラリが対応しているためである。下のtrain_img, test_imgはデータ数×28x28x1の4次元配列となっている。最後の「1」は色の数を表している(例えばRGBの3色の画像を学習させる場合はここに3が入る)。

In [38]:
train_img=x_train2.reshape(x_train2.shape[0],pixel_size,pixel_size,1)
test_img=x_test2.reshape(x_test2.shape[0],pixel_size,pixel_size,1)

これで、train_imgとtest_imgに28x28の画像ファイルが入る。次にパラメータを設定する。

In [39]:
num_classes=10
batch_size=128
epochs=30

num_classはカテゴリの数(MNISTの場合は0~9の10カテゴリ)である。batch_sizeは、全データを128個のカタマリに分けて CNNでは、全データを幾つかのサブセットに分けて学習するが、この幾つかに分けたぞれぞれのサブセットに含まれるデータの数をbatch_sizeと呼び、多くの場合に2のべき乗となっている。細かく分けたサブセットで学習を行い、それが1周して全データを用いた学習となるサイクルを1エポックと呼ぶ。学習をどれくらい繰り返すかを決める変数がepochsであり、今は30としている。

次に、Kerasでの学習に用いるcategoryの変数を用意する。この変数は、カテゴリ数×データ数の2次元配列で、それぞれのデータの正解カテゴリに対応する要素が1、そうでない要素が0の値が入っている。

In [40]:
y_train_category = keras.utils.to_categorical(y_train2, num_classes)
y_test_category = keras.utils.to_categorical(y_test2, num_classes)

下のセルのinput_shapeは入力画像のサイズを指定している。

さらに、model=Sequential()以降がCNNのモデルの構築の部分となっている。

model.addの中のConv2Dは畳み込みレイヤーを加えていることに対応している。MaxPooing2D, Flatten, Dropoutも同様にレイヤーを加えており、最後にmodel.compliで最適化に用いる関数などを指定することにより、モデルが構築される。

In [41]:
input_shape=(pixel_size, pixel_size,1)

model = Sequential()
model.add(Conv2D(32, kernel_size=(3, 3),
                     activation='relu',
                     input_shape=input_shape))
model.add(Conv2D(64, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='softmax'))

model.compile(loss=keras.losses.categorical_crossentropy,
                  optimizer=keras.optimizers.Adadelta(),
                  metrics=['accuracy'])
model.summary()
Model: "sequential_1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_1 (Conv2D)            (None, 26, 26, 32)        320       
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 24, 24, 64)        18496     
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 12, 12, 64)        0         
_________________________________________________________________
dropout_1 (Dropout)          (None, 12, 12, 64)        0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 9216)              0         
_________________________________________________________________
dense_4 (Dense)              (None, 128)               1179776   
_________________________________________________________________
dropout_2 (Dropout)          (None, 128)               0         
_________________________________________________________________
dense_5 (Dense)              (None, 10)                1290      
=================================================================
Total params: 1,199,882
Trainable params: 1,199,882
Non-trainable params: 0
_________________________________________________________________

この構築したモデルを実行するには、以下のようにする。

In [42]:
model.fit(train_img, y_train_category,
            batch_size=batch_size,
            epochs=epochs,
            verbose=1,
            validation_data=(test_img, y_test_category))
Train on 6000 samples, validate on 1000 samples
Epoch 1/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.9518 - accuracy: 0.6915 - val_loss: 0.3444 - val_accuracy: 0.8890
Epoch 2/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.3161 - accuracy: 0.9053 - val_loss: 0.1906 - val_accuracy: 0.9410
Epoch 3/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.2123 - accuracy: 0.9382 - val_loss: 0.1430 - val_accuracy: 0.9460
Epoch 4/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.1574 - accuracy: 0.9520 - val_loss: 0.1111 - val_accuracy: 0.9580
Epoch 5/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.1288 - accuracy: 0.9612 - val_loss: 0.1205 - val_accuracy: 0.9550
Epoch 6/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.1103 - accuracy: 0.9697 - val_loss: 0.0927 - val_accuracy: 0.9670
Epoch 7/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0902 - accuracy: 0.9725 - val_loss: 0.0870 - val_accuracy: 0.9710
Epoch 8/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0764 - accuracy: 0.9763 - val_loss: 0.0792 - val_accuracy: 0.9710
Epoch 9/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0680 - accuracy: 0.9793 - val_loss: 0.0734 - val_accuracy: 0.9750
Epoch 10/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0593 - accuracy: 0.9803 - val_loss: 0.0766 - val_accuracy: 0.9730
Epoch 11/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0487 - accuracy: 0.9850 - val_loss: 0.0728 - val_accuracy: 0.9740
Epoch 12/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0495 - accuracy: 0.9860 - val_loss: 0.0687 - val_accuracy: 0.9750
Epoch 13/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0329 - accuracy: 0.9887 - val_loss: 0.0768 - val_accuracy: 0.9690
Epoch 14/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0381 - accuracy: 0.9872 - val_loss: 0.0642 - val_accuracy: 0.9780
Epoch 15/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0334 - accuracy: 0.9885 - val_loss: 0.0587 - val_accuracy: 0.9770
Epoch 16/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0280 - accuracy: 0.9903 - val_loss: 0.0638 - val_accuracy: 0.9760
Epoch 17/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0266 - accuracy: 0.9910 - val_loss: 0.0649 - val_accuracy: 0.9780
Epoch 18/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0263 - accuracy: 0.9907 - val_loss: 0.0644 - val_accuracy: 0.9810
Epoch 19/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0240 - accuracy: 0.9923 - val_loss: 0.0651 - val_accuracy: 0.9780
Epoch 20/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0201 - accuracy: 0.9930 - val_loss: 0.0634 - val_accuracy: 0.9790
Epoch 21/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0217 - accuracy: 0.9927 - val_loss: 0.0558 - val_accuracy: 0.9820
Epoch 22/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0192 - accuracy: 0.9940 - val_loss: 0.0573 - val_accuracy: 0.9780
Epoch 23/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0196 - accuracy: 0.9938 - val_loss: 0.0547 - val_accuracy: 0.9830
Epoch 24/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0155 - accuracy: 0.9948 - val_loss: 0.0503 - val_accuracy: 0.9850
Epoch 25/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0144 - accuracy: 0.9955 - val_loss: 0.0576 - val_accuracy: 0.9810
Epoch 26/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0150 - accuracy: 0.9952 - val_loss: 0.0630 - val_accuracy: 0.9810
Epoch 27/30
6000/6000 [==============================] - 6s 1ms/step - loss: 0.0156 - accuracy: 0.9945 - val_loss: 0.0828 - val_accuracy: 0.9780
Epoch 28/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0158 - accuracy: 0.9950 - val_loss: 0.0585 - val_accuracy: 0.9840
Epoch 29/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0150 - accuracy: 0.9953 - val_loss: 0.0564 - val_accuracy: 0.9820
Epoch 30/30
6000/6000 [==============================] - 7s 1ms/step - loss: 0.0143 - accuracy: 0.9952 - val_loss: 0.0670 - val_accuracy: 0.9810
Out[42]:
<keras.callbacks.callbacks.History at 0x1c441c44ec8>

SVMやRFのときと同じにように、分類のaccuracyをチェックするためには、以下のようにする。2つの値が出力されるが、最初の値がtest dataのlossの値で、2番目の値がaccuracyである。

In [43]:
model.evaluate(test_img, y_test_category, verbose=0)
Out[43]:
[0.06698677124537061, 0.9810000061988831]

以下のようにすると、confusion matrixが表示されるはずである。

In [44]:
prediction = model.predict_classes(test_img)
classifier = confusion_matrix(y_test2, prediction)
plot_confusion_matrix(classifier,class_names)

発展問題 11-1:
上のCNNのモデル構築の部分で、model = Sequential()の行以降のmodel.add()の部分がCNNのレイヤーを追加している箇所である。今回用いたモデルは比較的単純なCNNであるが、より層を深くするなどによって、accuracyを上げることが出来るであろうか?例えばgoogleで「CNN MNIST チューン」とかを検索すると、いろいろなCNNが試されている。そうしたものを参考にして、MNISTやFashion MNIST、植物葉のデータなどの予測精度がモデルに応じてどのように変化するかを探索してみよ。