雑記

データを分析するとき、取り込んだままのデータの列をそのままつかうのではなくて、元々あった列をもとに新しい列をつくることとかあると思います。でも Python とかよくわからないので Google 検索してどうやって加工するか調べると思います。

例えば、元々あった列のうち a, b, c 列の中央値をとり、新しい列 d に追加したいと思います。
以下の (1) (2) はどちらも望みを満たすようであることがわかります。でも (1) は遅いのでやってはいけないとわかります。
というか (2) のように一気にできるとわかれば (1) のように行ごとの処理にはしないと思いますが調べ方が難しいと思います。

# -*- coding: utf-8 -*-
import pandas
import numpy

if __name__ == "__main__":
  data = pandas.DataFrame({'name': numpy.repeat('Alice', 100000),
                           'a': range(0, 100000, 1),
                           'b': range(0, 200000, 2),
                           'c': range(0, 400000, 4)})

  # (1) 44.674 sec
  data['d'] = data.apply(lambda x: numpy.median(x[['a', 'b', 'c']]), axis=1) 
  # (2) 0.010 sec
  data['d'] = numpy.median(data[['a', 'b', 'c']].values, axis=1)

ともあれ、numpy.mean や numpy.median などで各行への処理を一気にできるのだったらそうすればいいことはわかります。
ただ、もっと自分で定義した処理をしたいときはベクトル演算が定義されているとは限らないので困ると思います。
そのような場合は numpy.vectorize で関数をベクトル化することにより処理速度を改善できると以下の記事にあります。
sinhrks.hatenablog.com

例えば、データに元々ある time 列をもとに、米国サマータイム期間かどうかを is_dst 列として付加したいとします。
apply をつかうと (1) のようにかけます。でも上の記事のアドバイス通りこの処理をベクトル化した (2) の方が速いです。

# -*- coding: utf-8 -*-
import pandas
import numpy

def is_dst(dt):
  return(pandas.Timestamp(dt).tz_localize('Asia/Tokyo').tz_convert('US/Eastern').dst().seconds > 0)

if __name__ == "__main__":
  data = pandas.DataFrame({
    'name': numpy.repeat('Alice', 100000),
    'time': numpy.repeat(['2017/03/12 15:59', '2017/03/12 16:00'], [50000, 50000])})
  data.time = pandas.to_datetime(data.time)

  # (1) 28.834 sec
  data['is_dst'] = data.apply(lambda x:
    x.time.tz_localize('Asia/Tokyo').tz_convert('US/Eastern').dst().seconds > 0, axis=1)
  # (2) 9.768 sec
  data['is_dst'] = numpy.vectorize(is_dst)(data.time)

ところで、上の実装で自作の is_dst 関数は受け取った日付時刻を pandas.Timestamp 型に変換していますが、time 列は元々 pandas.Timestamp 型です。なのになぜこの変換をしているのかというと、ベクトル化した関数を time 列に適用したときに time 列が勝手に numpy.datetime64 型になってしまったからです。なんでそうなるのかは Python とかよくわからないのでわかりません。

雑記: sklearn.cluster.KMeans で k-means 法

sklearn.cluster.KMeans をつかってみたいときとかあると思います。iris データにつかうと以下のようになります。

# -*- coding: utf-8 -*-
import pandas
from sklearn import datasets
from sklearn.cluster import KMeans
from sklearn.metrics import confusion_matrix

if __name__ == "__main__":
  iris = datasets.load_iris()
  predict = KMeans(n_clusters=3).fit(iris.data)

  mat = pandas.DataFrame(confusion_matrix(iris.target, predict.labels_))
  mat.index = iris.target_names
  print(mat)
             0   1   2
setosa       0  50   0
versicolor  48   0   2
virginica   14   0  36

setosa はすべてクラスタ1に分類され(かつ、クラスタ1に他の品種は混ざっていない)、一方、virginica はクラスタ0とクラスタ2に分裂してしまったことがわかります。各品種が分類された最大クラスタが対角線上に並んでいないと見にくいという人は、以下のように列を入れ替えてください。

  mat = mat.loc[:,mat.idxmax(axis=1)]
  mat.columns = range(iris.target_names.shape[0])
  print(mat)
             0   1   2
setosa      50   0   0
versicolor   0  48   2
virginica    0  14  36

クラスタリングした結果、各クラスタの重心はどこになったかも確認できます(以下)。

  print(predict.cluster_centers_)
[[ 6.85        3.07368421  5.74210526  2.07105263]
 [ 5.006       3.418       1.464       0.244     ]
 [ 5.9016129   2.7483871   4.39354839  1.43387097]]

初期クラスタ重心を与えることもできます。
以下のようにすれば元データの 0 番目、50 番目、100 番目を初期クラスタ重心にできます。

# -*- coding: utf-8 -*-
import numpy
import pandas
from sklearn import datasets
from sklearn.cluster import KMeans
from sklearn.metrics import confusion_matrix

if __name__ == "__main__":
  iris = datasets.load_iris()
  initial_centers = numpy.array([iris.data[0], iris.data[50], iris.data[100]])
  predict = KMeans(n_clusters=3, init=initial_centers).fit(iris.data)

  mat = pandas.DataFrame(confusion_matrix(iris.target, predict.labels_))
  mat.index = iris.target_names
  print(mat)
  print(predict.n_iter_)

ただし、これを実行すると以下のように警告文が出てきます。これの意味は、KMeans はクラスタ初期化方法が k-means++(デフォルト)や random の場合はデフォルトで n_init=10 通りの初期値からのクラスタリングを試み、最もよいクラスタリング結果を出力します(モーメントがよく更新され続けた回をもっともよいクラスタリング結果とする)。しかし、初期クラスタ重心が明示的に与えられたときにはその他の初期化を試みるのはおかしいので、デフォルトの n_init=10 は無視して1回しかクラスタリングを試みませんよという意味です。その代わりにというか、クラスタリングイテレーション回数 n_iter_ が返却されます。クラスタ重心は3回更新されるのみだったようです。

/usr/local/lib/python2.7/site-packages/sklearn/cluster/k_means_.py:889: RuntimeWarning:
 Explicit initial center position passed: performing only one init in k-means
 instead of n_init=10
  return_n_iter=True)
             0   1   2
setosa      50   0   0
versicolor   0  48   2
virginica    0  14  36
3

雑記: Keras で自作レイヤー(Dense)

Keras でオリジナルの自作レイヤーを追加したいときとかあると思います。
自作レイヤー自体は以下の記事でつかったことがありますが、これはウェイトをもつレイヤーではなく、最後にかぶせて損失関数のみをカスタマイズするためのレイヤーでした。
Keras で変分自己符号化器(VAE)を学習したい - クッキーの日記
ウェイトをもつレイヤーはどう追加するのか知りたいのですが、以下のドキュメントの通りにすればできるのか確認します。
Writing your own Keras layers - Keras Documentation
上のページの MyLayer は単純な Dense のようなので、実際に Dense に置き換えられるか確認します。

keras.layers.Dense をつかった場合

まず、 普通に Keras の Dense をつかって MNIST の分類器をつくると例えば以下になると思います。モデルは適当です。
→ 下記スクリプトの設定で 95% 程度のテスト精度になりました。

# -*- coding: utf-8 -*-
import numpy
from keras.layers import Input, Dense, Dropout
from keras.models import Model
from keras.datasets import mnist
from keras.utils import np_utils

if __name__ == "__main__":
  (x_train, y_train), (x_test, y_test) = mnist.load_data()
  x_train = x_train.astype('float32') / 255.
  x_test = x_test.astype('float32') / 255.
  x_train = x_train.reshape((len(x_train), numpy.prod(x_train.shape[1:])))
  x_test = x_test.reshape((len(x_test), numpy.prod(x_test.shape[1:])))
  y_train = np_utils.to_categorical(y_train, num_classes=10)
  y_test = np_utils.to_categorical(y_test, num_classes=10)

  input = Input(shape=(x_train.shape[1],))
  converted = Dense(20, activation='relu')(input)
  converted = Dropout(0.2)(converted)
  converted = Dense(20, activation='relu')(converted)
  converted = Dropout(0.2)(converted)
  converted = Dense(10, activation='softmax')(converted)
  classifier = Model(input, converted)
  classifier.compile(optimizer='adam', loss='categorical_crossentropy',
                     metrics=['accuracy'])
  print(classifier.summary())

  classifier.fit(x_train, y_train, epochs=20, batch_size=100, shuffle=True)

  scores = classifier.evaluate(x_train, y_train, verbose=0)
  print('Accuracy (train): %.2f%%' % (scores[1]*100))
  scores = classifier.evaluate(x_test, y_test, verbose=0)
  print('Accuracy (test) : %.2f%%' % (scores[1]*100))
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         (None, 784)               0         
_________________________________________________________________
dense_1 (Dense)              (None, 20)                15700     
_________________________________________________________________
dropout_1 (Dropout)          (None, 20)                0         
_________________________________________________________________
dense_2 (Dense)              (None, 20)                420       
_________________________________________________________________
dropout_2 (Dropout)          (None, 20)                0         
_________________________________________________________________
dense_3 (Dense)              (None, 10)                210       
=================================================================
Total params: 16,330
Trainable params: 16,330
Non-trainable params: 0
_________________________________________________________________
Accuracy (train): 95.69%
Accuracy (test) : 94.92%
自作 Dense をつかった場合

次に、import Dense をつかわずに自作 Dense ( class MyDense ) に置き換えたのが以下です。
→ これも 95% 程度のテスト精度になったので、ちゃんと Dense になっていると思います。
  パラメータ数も上と同じ 16,330 です。

Keras ドキュメント(Writing your own Keras layers - Keras Documentation)にあるとおり、自作レイヤーのクラスでは build メソッドでウェイトを定義し、call メソッドにウェイトを用いたこのレイヤーの計算処理を実装すればよいはずです。compute_output_shape メソッドは次のレイヤーのために出力 shape を返すようにしておきます。

# -*- coding: utf-8 -*-
import numpy
from keras import backend as K
from keras.engine.topology import Layer
from keras.layers import Input, Dropout
from keras.models import Model
from keras.datasets import mnist
from keras.utils import np_utils

class MyDense(Layer):
  def __init__(self, output_dim, activation, **kwargs):
    self.output_dim = output_dim
    self.activation = activation
    super(MyDense, self).__init__(**kwargs)

  def build(self, input_shape):
    self.kernel = self.add_weight(name='kernel',
                                  shape=(input_shape[1], self.output_dim),
                                  initializer='glorot_uniform')
    self.bias   = self.add_weight(name='bias',
                                  shape=(1, self.output_dim),
                                  initializer='zeros')
    super(MyDense, self).build(input_shape)

  def call(self, x):
    if self.activation == 'relu':
      return(K.relu(K.dot(x, self.kernel) + self.bias))
    elif self.activation == 'softmax':
      return(K.softmax(K.dot(x, self.kernel) + self.bias))
    else:
      return(K.dot(x, self.kernel) + self.bias)

  def compute_output_shape(self, input_shape):
    return(input_shape[0], self.output_dim)

if __name__ == "__main__":
  (x_train, y_train), (x_test, y_test) = mnist.load_data()
  x_train = x_train.astype('float32') / 255.
  x_test = x_test.astype('float32') / 255.
  x_train = x_train.reshape((len(x_train), numpy.prod(x_train.shape[1:])))
  x_test = x_test.reshape((len(x_test), numpy.prod(x_test.shape[1:])))
  y_train = np_utils.to_categorical(y_train, num_classes=10)
  y_test = np_utils.to_categorical(y_test, num_classes=10)

  input = Input(shape=(x_train.shape[1],))
  converted = MyDense(20, activation='relu')(input)
  converted = Dropout(0.2)(converted)
  converted = MyDense(20, activation='relu')(converted)
  converted = Dropout(0.2)(converted)
  converted = MyDense(10, activation='softmax')(converted)
  classifier = Model(input, converted)
  classifier.compile(optimizer='adam', loss='categorical_crossentropy',
                     metrics=['accuracy'])
  print(classifier.summary())

  classifier.fit(x_train, y_train, epochs=20, batch_size=100, shuffle=True)

  scores = classifier.evaluate(x_train, y_train, verbose=0)
  print('Accuracy (train): %.2f%%' % (scores[1]*100))
  scores = classifier.evaluate(x_test, y_test, verbose=0)
  print('Accuracy (test) : %.2f%%' % (scores[1]*100))
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         (None, 784)               0         
_________________________________________________________________
my_dense_1 (MyDense)         (None, 20)                15700     
_________________________________________________________________
dropout_1 (Dropout)          (None, 20)                0         
_________________________________________________________________
my_dense_2 (MyDense)         (None, 20)                420       
_________________________________________________________________
dropout_2 (Dropout)          (None, 20)                0         
_________________________________________________________________
my_dense_3 (MyDense)         (None, 10)                210       
=================================================================
Total params: 16,330
Trainable params: 16,330
Non-trainable params: 0
_________________________________________________________________
Accuracy (train): 95.83%
Accuracy (test) : 95.14%

class MyDense の実装は Keras ドキュメント(Writing your own Keras layers - Keras Documentation)の MyLayer を参考にしましたが、必要に応じて keras.engine.topology.Layer のソースを参照して修正しました。
keras/topology.py at master · fchollet/keras · GitHub

  • 当初 Keras ドキュメントの MyLayer のとおりに実装したところ、パラメータ数からバイアス項がないことに気付いたので、バイアス項を追加しました。
  • add_weight の引数に name が必要だったので追加しました。name をつけて何につかうのかわかりません。→ 後から自作 MyDense にはバイアス項がないことに気付いたとき、ユニークな名前でいくつでも weight 行列を追加できることに気付きました。
  • kernel と bias の初期化方法は本家 Dense に合わせました(初期化が uniform だといまいち到達精度が 94% 程度にしかなりませんでした)(誤差もあるかもしれませんが)。
  • 今回のモデルでは Dense を ReLU または softmax で活性化するので、Dense よろしく MyDense インスタンス生成時に活性化関数も指定するようにしました(ReLU と softmax 以外未対応)(活性化関数自体を引数に取ればいいんだろうとは思います)。



これで Dense が自作できることは確認できたのですが、自分がつくりたいのは Dense ではないのでこれからもっと調べていきます。

雑記: 対象日付が米国サマータイム期間か判定(R, Python)

時刻カラムがあるようなデータのプレ処理で、各レコードが米国サマータイム期間かどうか自動判定したいときとかあると思います。サマータイムかそうでないかで東京とニューヨークの時差が違うからね。
とりわけ、時刻カラムのタイムゾーンは日本時間なのですが、その時刻がニューヨークでサマータイムかどうかがほしいことがあると思います。以下の IS_DST 列を追加したいということです(サマータイム突入時には、現地では深夜2時を3時にジャンプさせますが、このジャンプ時刻は日本時間だと16時になります)。
DATETIMEIS_DST
2017/03/12 15:59:590
2017/03/12 16:00:001
時刻文字列を入れると、日本時間のその時刻にニューヨークではサマータイムかどうか返す関数があればよく、R だと例えば以下のようにかけます(参考: R の日付時刻オブジェクト - クッキーの日記)。

is_dst <- function(datetime_string='2017/03/12 15:59:59') {
  time_lt <- strptime(datetime_string, format='%Y/%m/%d %H:%M:%S')
  time_ct <- as.POSIXct(time_lt, tz='Japan')
  return(as.POSIXlt(time_ct, tz="EST5EDT")$isdst == 1)
}

print(is_dst('2017/03/12 15:59:59')) # FALSE
print(is_dst('2017/03/12 16:00:00')) # TRUE

ミソは、isdst のようなプロパティを取り出せるのは POSIXlt 型ですが、POSIXlt 型のままタイムゾーン間を行き来することはできないので、一旦 POSIXct 型を経由するところです(もっとスマートな方法があるかもしれないですがわからない)。

一方、Python では例えば以下のようにかけます(だいたい pandas をつかうので pandas.Timestamp 型をつかいます)。

# -*- coding: utf-8 -*-
import pandas

def is_dst(datetime_string='2017/03/12 15:59:59'):
  datetime = pandas.Timestamp(datetime_string)
  datetime = datetime.tz_localize('Asia/Tokyo')
  datetime = datetime.tz_convert('US/Eastern')
  return(datetime.dst().seconds > 0)

if __name__ == "__main__":
  print(is_dst('2017/03/12 15:59:59')) # False
  print(is_dst('2017/03/12 16:00:00')) # True

時刻文字列をタイムスタンプに変換して、タイムゾーンを日本時間にして、アメリカ東部標準時に変換するだけです。dst() というメソッドは Standard Time との時差(つまり、サマータイムのときは1時間で、そうでなければ0時間)が取得できるようなので、0 より大きければサマータイムです。



上の Python スクリプトでいま家では上手くいったんだけど、会社ではある1つの日付文字列から上の方法でサマータイムかどうかを上手く取得できなかった気がするんだよね(何かエラーが出た)。Pandas DataFrameで日時データのタイムゾーン変換 - pLog という記事によると、データフレーム型ではタイムスタンプをインデックスにしないと tz_convert ができないとあります。となると、ある1つの文字列に対して判定を行いたい場合、以下のようにすべきなのでしょうか…?

def is_dst(datetime_string='2017/03/12 15:59:59'):
  df = pandas.DataFrame({'datetime': [datetime_string]})
  df.index = pandas.DatetimeIndex(df.datetime, name='datetime')
  df.index = df.index.tz_localize('Asia/Tokyo')
  df.index = df.index.tz_convert('US/Eastern')
  return(df.index[0].dst().seconds > 0)

というかある1つの文字列に対して判定を行いたいだけなら pandas をつかわなくてもいいですね。
Python: How do you convert datetime/timestamp from one timezone to another timezone? - Stack Overflow

# -*- coding: utf-8 -*-
import datetime
import pytz

def is_dst(datetime_string='2017/03/12 15:59:59'):
  dt0 = datetime.datetime.strptime(datetime_string, '%Y/%m/%d %H:%M:%S')
  tz0 = pytz.timezone('Asia/Tokyo')
  tz1 = pytz.timezone('US/Eastern')
  dt1 = tz0.localize(dt0).astimezone(tz1)
  return(dt1.dst().seconds > 0)

Keras で変分自己符号化器(VAE)を学習したい

以下の記事の続きです。Kerasブログの自己符号化器チュートリアルをやるだけです。
Keras で自己符号化器を学習したい - クッキーの日記



Kerasブログの自己符号化器チュートリアルBuilding Autoencoders in Keras)の最後、Variational autoencoder(変分自己符号化器;VAE)をやります。VAE についてのチュートリアル上の説明は簡単なものなので、以下では自分で言葉を補っています。そのため、不正確な記述があるかもしれません。

変分自己符号化器(VAE)って何

そのデータが生成するメカニズムに仮定をおいているとき(そのデータの生成モデルを仮定しているとき)、モデルのパラメータの最適化をするのに VAE を用いることができます。今回は、「それぞれの手書き数字には、その手書き数字に対応する隠れ変数の確率分布があって、その分布からの実現値を変換することによって手書き数字ができている」と仮定しています。さらに、その隠れ変数の次元は2次元であって(  \in {\mathbb R}^2 )、隠れ変数の確率分布は2次元正規分布 \mu, \, \Sigma によってのみ決まる )と仮定しています。この分布からの実現値が f: {\mathbb R}^2 \to {\mathbb R}^{784} なる写像によって 784 ピクセルに変換されて目の前に手書き数字として具現化したのだと考えているわけです。

f:id:cookie-box:20170513202623p:plain:w360
なので、今回訓練するのは以下のようなモデルになります。
後から気付いたのですが、分散ベクトルの次元が2であることからわかるように、第1変数と第2変数が独立な2次元正規分布を考えているようです。したがって、図のように分布は斜めに傾きません。が、面倒なので描き直しません。

f:id:cookie-box:20170513200541p:plain:w720

ここまでの自己符号化器では、入力と出力の交差エントロピーを損失関数としてモデルを訓練していました。VAE でも同様に訓練してもよいのですが、今回は入力と出力の交差エントロピーに、隠れ変数空間の事前分布と事後分布のKL情報量も加えたものを損失関数とします。つまり、隠れ変数の確率分布が大幅に更新されることに制約を課すわけですが、こうすることで、隠れ変数をきれいに( well-formed ? )モデリングでき、過学習も防げるそうです。

  • このような特殊な損失関数を実現するために、訓練時には自分で定義したレイヤーを最後にかぶせて、そのレイヤーで入出力の交差エントロピーと事前事後分布のKL情報量を損失関数に加えています(スクリプトは記事の一番最後)。そして、compile() では損失関数を指定していません( loss=None )。Keras の仕様をきちんと確認していないのですが、Layerクラスを継承して自分で定義したレイヤーでは損失関数を付加するということができ、それより手前の層に誤差逆伝播するのだと思います。たぶん。
実行結果

今回、上図の VAE を MNIST の 60000 の訓練用データで学習しました(バッチサイズ:100、エポック数:10)。
テストデータについて、隠れ変数の確率分布の平均をプロットすると以下のようになりました。同じ数字のデータを同じ色でプロットしています。左上の濃い紫のかたまりが "0" です。左下のもう少し薄い紫のかたまりが "1" です。右下の黄緑色のかたまりは "7" です。同じ数字が隠れ変数空間上でかたまる傾向が確認できます。ただ、Kerasブログの元記事とは数字のかたまり方の形や配置が結構異なっています。

   f:id:cookie-box:20170513223659p:plain:w550

また、隠れ変数空間から均等にサンプリングした点がどんな「手書き数字」を生成するか(どうデコードされるか)というのも確かめることができます。以下のようになりました。実行したスクリプトの関係で上のプロットと上下が反転していますが、上のプロットと対応するように左上の方が "1"、右上の方が "7" に見えるのが確認できます。下の方は "0" に見えます。右側の真ん中あたりが "9" に見えますが、ちょうど上の図でも右側の真ん中あたりに "9" (黄色)が見えます。
f:id:cookie-box:20170513223749p:plain:w770

スクリプト

補助関数の定義とクラスの定義を冒頭にもってきているのと日本語のコメントを加えている以外は以下にあるスクリプトと同じです。ただし、2番目のプロットはブログ記事のスクリプトと記述が違ったのでブログ記事の方に合わせています。
keras/variational_autoencoder.py at master · fchollet/keras · GitHub

# -*- coding: utf-8 -*-
import numpy
import matplotlib.pyplot as plt
from scipy.stats import norm
import sys

from keras.layers import Input, Dense, Lambda, Layer
from keras.models import Model
from keras import backend, metrics
from keras.datasets import mnist

# 2次元正規分布から1点サンプリングする補助関数です
def sampling(args):
  z_mean, z_log_var = args
  epsilon = backend.random_normal(shape=(batch_size, latent_dim), mean=0., stddev=epsilon_std)
  return z_mean + backend.exp(z_log_var / 2) * epsilon

# Keras の Layer クラスを継承してオリジナルの損失関数を付加するレイヤーをつくります
class CustomVariationalLayer(Layer):
  def __init__(self, **kwargs):
    self.is_placeholder = True
    super(CustomVariationalLayer, self).__init__(**kwargs)

  def vae_loss(self, x, x_decoded_mean): # オリジナルの損失関数
    # 入力と出力の交差エントロピー
    xent_loss = original_dim * metrics.binary_crossentropy(x, x_decoded_mean) 
    # 事前分布と事後分布のKL情報量
    kl_loss = - 0.5 * backend.sum(1 + z_log_var - backend.square(z_mean) - backend.exp(z_log_var), axis=-1)
    return backend.mean(xent_loss + kl_loss)

  def call(self, inputs):
    x = inputs[0]
    x_decoded_mean = inputs[1]
    loss = self.vae_loss(x, x_decoded_mean)
    self.add_loss(loss, inputs=inputs) # オリジナルの損失関数を付加
    return x # この自作レイヤーの出力を一応定義しておきますが、今回この出力は全く使いません

if __name__ == "__main__":
  batch_size = 100
  original_dim = 784
  latent_dim = 2
  intermediate_dim = 256
  epochs = 10
  epsilon_std = 1.0

  #============================================================
  # 変分自己符号化器を構築します

  # エンコーダ
  x = Input(batch_shape=(batch_size, original_dim))
  h = Dense(intermediate_dim, activation='relu')(x)
  z_mean = Dense(latent_dim)(h)
  z_log_var = Dense(latent_dim)(h)
  z = Lambda(sampling, output_shape=(latent_dim,))([z_mean, z_log_var])

  # デコーダ
  decoder_h = Dense(intermediate_dim, activation='relu')
  decoder_mean = Dense(original_dim, activation='sigmoid')
  h_decoded = decoder_h(z)
  x_decoded_mean = decoder_mean(h_decoded)

  # カスタマイズした損失関数を付加する訓練用レイヤー
  y = CustomVariationalLayer()([x, x_decoded_mean])
  vae = Model(x, y)
  vae.compile(optimizer='rmsprop', loss=None)

  #============================================================
  # モデルを訓練します

  (x_train, y_train), (x_test, y_test) = mnist.load_data()
  x_train = x_train.astype('float32') / 255.
  x_test = x_test.astype('float32') / 255.
  x_train = x_train.reshape((len(x_train), numpy.prod(x_train.shape[1:])))
  x_test = x_test.reshape((len(x_test), numpy.prod(x_test.shape[1:])))
  vae.fit(x_train, shuffle=True, epochs=epochs, batch_size=batch_size,
          validation_data=(x_test, x_test))

  #============================================================
  # 結果を表示します

  # (1) 隠れ変数空間のプロット(エンコードした状態のプロット)
  encoder = Model(x, z_mean) # エンコーダのみ分離
  x_test_encoded = encoder.predict(x_test, batch_size=batch_size)
  plt.figure(figsize=(6, 6))
  plt.scatter(x_test_encoded[:, 0], x_test_encoded[:, 1], c=y_test)
  plt.colorbar()
  plt.show()

  # (2) 隠れ変数空間からサンプリングした点がどんな手書き数字を生成するか(どうデコードされるか)をプロット
  decoder_input = Input(shape=(latent_dim,))
  _h_decoded = decoder_h(decoder_input)
  _x_decoded_mean = decoder_mean(_h_decoded)
  generator = Model(decoder_input, _x_decoded_mean) # デコーダのみ分離

  n = 15  # 15x15 個の手書き数字をプロットする
  digit_size = 28
  figure = numpy.zeros((digit_size * n, digit_size * n))
  grid_x = numpy.linspace(-15, 15, n)
  grid_y = numpy.linspace(-15, 15, n)

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

  plt.figure(figsize=(10, 10))
  plt.imshow(figure)
  plt.show()