故障予測問題に TensorFlow Probability を利用する(モデル推定編)

15 min readJan 6, 2019

この記事は TensorFlow Advent Calendar 2018 の6日目にすべりこませるために書いたものですが、先に枠が埋まってしまったので冬休みの宿題ということにしました。5日目にすべりこませた「データ生成編」の続きです。

TensorFlow Probability は TensorFlow をベースにした統計モデルを構築するための Python ライブラリです。モデル構築に用いる要素として確率変数、分散関数(distributions)、bijectors という確率分布を変形する仕組みを持ち、モデル推定のために MCMC や変分推定(Variational Inference)の機能を持ちます。またいわゆる Bayesian Neural Network を構築するための部品も備えています。

筆者は Bayesian の考えかたにも TensorFlow Probability (以下TFP)の利用法にも不慣れなのですが、Bayesian Methods for Hackers というプログラマ向けに Bayesian をわかりやすく解説した書籍の TFP 版が公開されたのをきっかけにして両方に一度に入門してみました。

長くなりそうなので「データ生成編」と「モデル推定編」に分けました。今回は「モデル推定編」です。

追記(2019–01–10)

記事公開後 HELLO CYBERNETICS さんにTwitterでコメントをいただき、このコードで行なっているのは変分推論によるNNのパラメーター(kernel, weights)の事後分布推定であるとわかりました。

また「点推定」という語彙は通常モデルのパラメーター(よくθと表記される。今回のケースだと kernel, weigts 自体)自体を確率変数ではなく具体的な値に決めるということを指す場合に使われるようでしたので、別の単語を使うように表現を修正しました。

モヤっとしていたところがこのやりとりのなかですっきりして、非常にありがたいコメントでした。ありがとうございます。

問題設定(再び)

問題設定として工場の装置の故障予測という仮想のタスクを考えます。元となるアイデアは “A strategy for implementing industrial predictive maintenance” という一連のブログ記事から得ています。工場の装置のメンテナンス計画の最適化のために装置の故障予測を行いたい、というシナリオです。このブログ記事内ではタスクとしてある装置が一定期間内に故障するか否かの判定(分類)、寿命がどのくらい残っているかの予測(回帰)、装置の挙動の異常の検知(異常検知)、装置のパラメーターなどの操作の学習(強化学習)といった手法があるとしていますが、Bayesian Methods を適用するにあたり、最もシンプルな分類のタスクをベースに以下のように故障確率の推定というタスクに変更します。

ある装置が今後一定期間内に故障するベルヌーイ分布のパラメータpの確率分布を求める
故障する/しないの2値なのでベルヌーイ分布に従うとする
(つまりある値p(0 ≤ p ≤ 1)をパラメータとして確率pで故障する、(1-p) で故障しないとする)
上記のpの値の取る確率分布を求める

「データ生成編」に書いたように予測に使える入力パラメーターとは別に直接知ることができないパラメーターがあり、確率分布も確実には知ることができないという前提になっています。

このようなタスク設定にすることを思いついたのは、TensorFlow Probability のフォーラムのこちらの投稿を読んだためです。この投稿で述べられているaleatoric uncertaintyとepistemic uncertaintyという概念を自分なりに解釈して上のように2段階の不確実性を持つモデル化をすることで、後に述べるような方法で生成するデータにも対応できるのではないかと考えました。(あくまでこのフォーラム投稿に着想を得たというだけで aleatoric uncertaintyとepistemic uncertainty という概念が本当にここで述べたような意味なのか自信ありません)

さて「データ生成編」で用意したのは配管に使うパイプの連結器具の、以下の入力パラメーター毎の故障の有無(0 or 1)を含む CSV ファイルでした。

稼動日数(延べ何日利用したか)
最後にメンテナンスしてからの日数
パイプ内の温度センサー情報

Bayesian Methods for Hackers (の最初のほう)ではデータを分析してPoisson分布とか名前のついている分布でモデル化して、その分布のパラメーターの確率分布をMCMCを用いて推定する、ということをしているのでまずはこのデータがどんな形をしているのかみてみます。

データ可視化

まずは pandas で CSV を読み込み基本的な統計情報をみてみます(詳細はNotebook参照)。故障のデータが全体の6.65%と少なめで偏りがあるのが注意が必要そうです。

次に稼動日数ごとに故障のデータの率をグラフにしてみます。

「データ生成編」の生成方法をみかえすとわかりますが、稼動日数が多いデータは少なめになってるので、日数多めのほうはイレギュラーの影響もありそうですが、おおむね日数が増えるに従って故障頻度が増えていることはわかります。

次は同様に温度の影響をグラフ化します。

これも温度が高くなるほど故障頻度が増えていることがみてとれます。

メンテナンス後日数についても同様にグラフにしてみましたが、これはあまり結果に主要な影響はなさそうで、強いていえば長いほうが頻度は上がる傾向がありそう、という程度でした。

そこで主要な影響を持ってそうな稼動日数と温度から故障頻度への3Dグラフを描画してみました。

うーん、ともかく

温度が高い(100を超えるあたりから)と稼動日数が0付近でも故障頻度は高い
温度が低くても稼動日数が長くなると故障頻度が高くなる

ということはなんとなくわかります。

モデル構築

さてここで3つの入力を元に故障確率pの分布を求めるモデルを Bayesian Methods for Hackers の例のように考えようと思っていたのですが、正直いって多変数の時にどのようなモデル化をしたらいいのかわかりませんでした。
また実践的にも対象のデータの事前知識をあまり必要としないで使える手法が欲しかったので、ここは Neural Network のweightsとbiasに確率変数を導入した Bayesian Neural Network を使ってみることにします。
Bayesian Neural Network の実装は tensorflow のリポジトリに含まれているサンプル([1], [2])を参考にしました。

net = tf.keras.Sequential([
    tfp.layers.DenseFlipout(units=20, activation=tf.nn.relu),
    tfp.layers.DenseFlipout(units=20, activation=tf.nn.relu),
    tfp.layers.DenseFlipout(units=1, activation=tf.nn.sigmoid),
])logits = net(features)
labels_distribution = tfp.distributions.Bernoulli(probs=logits)

モデル構築は tf.keras の API を利用して、レイヤーとして tfp.layers.DenseFlipout()を中間層に2層(活性化関数はReLU)、出力層としてユニット数1の層(活性化関数はsigmoid)を起きます。
出力層の出力 0<p<1 が Bernoulli 分布の p になります。
DenseFlipout は Dense と同様の全結合層ですが、weightsとbiasがただの変数ではなく確率変数である点が違います。同じ入力に対しても評価するたびに異なる p が出力されます。
なお DeneFliout 層は現状 eager execution に対応していないらしくてトレーニング時にエラーになってしまうので、graph mode を利用しています。

neg_log_likelihood = -tf.reduce_mean(
    labels_distribution.log_prob(tf.reshape(labels, (-1, 1))))
kl = sum(net.losses) / len(df_over)
loss = neg_log_likelihood + kl

loss 関数は通常のNNとは少し違っています。特に KL divergence の項が追加されるところが特徴的です。実装をちらっとみた限りだと通常の正規化項に似てるのかな……という感じですがよくわかってません。このあたりはサンプルの通り。

optimizer = tf.train.AdamOptimizer(learning_rate=0.00001)
train_op = optimizer.minimize(loss)init_op = tf.group(
        tf.global_variables_initializer(),
        tf.local_variables_initializer())with tf.Session() as sess:
  sess.run(init_op)
  for step in range(10001):
    _ = sess.run([train_op, accuracy_update_op])
    if step % 500 == 0:
    loss_value, accuracy_value = sess.run([loss, accuracy])
    print("Step: {:>3d} Loss: {:.3f} Accuracy: {:.3f}".format(step, loss_value, accuracy_value))
        
  tf.saved_model.simple_save(sess, export_dir="./model", inputs={"features": features}, outputs={"logits": logits, "labels_sample": labels_distribution.sample()})

ここがちょっと悩んだところで、BNNのトレーニングは勾配法のOptimizerを使って行ないます。
Bayesian Methods for Hackers の例だとモデルのパラメーター(BNNの場合はweightsとbiasの分布)は事前分布は指定する(Gaussian 分布)ものの事後分布にはなにか名前のついた特定の形状の分布を前提とせずに確率分布を推定するのにMCMCを利用しますが、DenseFlipoutの実装を確認してみたところ、weightsとbiasの事後分布はNormal分布に従うという前提になっていて(近似されていて)、そのパラメーターの確率分布の分布関数であの Normal の mean, scale は特定の値になります。Bayesian Methods for Hackers で例にあがっている手法だとこの場合 weights/bias といったパラメーターの分布をMCMCで推定するということになると思うのですが。
最初この違いがどういう意味をもつのか良くわからなかったのですが、冒頭に追記したようにコメントを頂いて、これが変分推論という手法による推論で、MCMCとは別の手法ではあるけどモデルのパラメーターの事後分布を推論しているという意味では同じと理解しました。ありがとうございました。

トレーニング結果

最終的な Accuracy は約 93.7% (実行するたびに変化するので目安ですが)と良くない印象です。Recall と F値も計算してみるとそれぞれ 0.12, 0.2 前後の値なので、やはり0(故障なし)の出力に偏っているようです。

得られたモデルは確率変数を出力するので、いくつか条件を固定してサンプリングして可視化してみます。

メンテナンス/温度を固定して、稼動日数を変化させた時に100回ずつサンプリングして散布図にしたものです。稼動日数が小さいうちはp<0.2の範囲に固まってますが、稼動日数が増えるに従って分散が大きくなっています。データ可視化のところにも書いたように生成したダミーデータは稼動日数が大きいデータ数が少ないので、その影響もあるかもしれません。

次は稼動日数/メンテナンスは固定して温度を変化させた時の散布図です。今度は温度が高くなるにつれて p が 0 付近に収束していくという分布になってしまって、これは期待と反します。温度についての学習がなぜうまくいかなかったのかはよくわかりません。

最後に稼動日数/メンテナンス/温度全て固定して、出力されるpの分布をみてみます。稼動日数140日の時の分布です。

2つのピークがある分布になりました。ここで「データ生成編」での神の視点を思い出すと、実はパイプには製造ラインという隠れたパラメーターがあり、稼動日数による故障確率のしきいちがずれています(140日と180日)。140日頃にはこの2つのグループの故障率は乖離していて、2つの山が見えることを期待していたのですが、実際に入力には陽に製造ラインの情報がないにもかかわらずこのピークを再現できたのにはちょっと驚きました。

はまったポイント

Overampling

データは故障時のものが少ないinbalancedなものだったので、故障(label=1)のものをoversamplingして比率を等しくしてみようとしたのですが、そうするとトレーニング時にかなりの頻度でlossがNaNになってしまう(おそらくweights/biasのGaussian分布のscaleが0になってしまうのではないかと思います)という現象がみられたのでやめました。

BNNの中間層ユニット数

最初は2層の中間層のユニット数を20にしていましたが、それだとlearning_rateを調整してもなぜかうまく収束せず、10にすると比較的にうまく(あくまで比較的ですが)いくようになりました。
またこれはただの印象ですがネットワークのサイズとlearning rateの関係の調整がシビアで通常のNNの時のカンが通じにくい気がしました。

tfp.distribution.Bernoulli() と log_prob() の引数のshape

tensorflow-probability での実装上でのpitfallですが tfp.distribution.Bernoulli() に渡す probs 引数と、その log_prob() に渡す引数の shape が一致していないと、意図しないbroadcastingが起きてloss計算が間違っていてうまく収束しないのになかなか気がつかなくてかなり時間を使ってしまいました。

まとめ

「データ生成編」で作ったダミーデータを元にして、TensorFlow Probability による Bayesian Neural Network を実装し学習しました。

最終的な予測性能はあまり良いとは言えませんでしたがモデルやトレーニングのパラメーターの調整不足なのかもしれません。通常のMLPを使った点推定の成績とも比較してみたいです。

またBNNの一般的なものなのか、TensorFlow Probability の現状の実装によるものか、わたしの実装がまずいからなのかはわかりませんが、モデルのサイズや learning rate などのパラメーター調節が難しいという印象でした。

予測性能には不満があるものの一方で、もうひとつの目的だった「データに隠された構造を発見できるか?」という疑問は、隠れパラメーターである製造ラインの影響による確率分布の2つのピークを再現できたので成功した、ような気がします。確率分布をみることで、それが何によるものかは不明ながら、パイプの故障率を左右するなにかの要因があるということが発見できるのではないかと思います。この解釈で正しいのかはちょっと自信がありませんが……。

実験に使ったノートブックは以下の Gist に保存しました。記事や実装の間違いなど指摘はコメントやTwitterアカウント(@nagachika) にお願いします。