はじめに
TensorFlow Probabilityには様々な確率分布が実装されています。 そして確率分布のパラメータをTensorFlowのVariable型で置いてやり、対数尤度関数を書き下してやれば、自動微分機能を駆使して簡単にパラメータ推定を実行できてしまいます。
今回はガウシアンプロセスの最尤推定を実行してみましょう。
ガウシアンプロセスで最低限知ってほしいこと
ガウシアンプロセスというのはまともに勉強すると、かなり難しいです。
しかし、とりあえずガウシアンプロセスという確率分布（確率過程）があるらしいと認めて、そのパラメータを推定すれば良いんだ！と思ってしまえば、ツールとして使う上ではそれほど難しいことを考える必要はありません。
線形回帰
ガウシアンプロセスで回帰をざっくり説明するため、まずは普通の線形回帰から見てみます。入力変数 $\bf x$、パラメータ$\bf w$として出力変数 $y$が下記で表されるものが普通の線形回帰です。
$$ \begin{align} y &= \bf {w^T x} + \epsilon \\ \epsilon &\sim \mathcal N (0, \sigma^2) \end{align} $$
ちなみに$\epsilon$は平均 $0 $ 分散 $\sigma^2 $ のガウス分布から生起しているノイズという扱いです。 \epsilonが確率変数でガウスノイズであるとすると、出力である確率変数$y$もガウス分布から生起しているとみなせます。具体的には下記のように書き下されます。
$$ y \sim \mathcal N (\bf {w^T x}, \sigma^2) $$
これは平均的には $\bf w^T x$ なのだがガウス分布から生起しているノイズ $\epsilon$ の影響で $y$ がガウス分布から生起しているように振る舞うという意味になります。こいつを最尤推定して $w$ を求めるのがいわゆる線形回帰であり、最小二乗法と呼ばれる手法です。
ガウシアンプロセス回帰
さて、ガウシアンプロセス回帰は普通の線形回帰に比べて何が違うのかというと、上記の線形回帰と比較すると意外と単純で
$$ y \sim \mathcal N (f({\bf x}), \sigma^2) $$ となります。$\bf {w^T x}$ が $f({\bf x})$ に代わりました！なんということでしょうか。
私たちはパラメータ $ \bf w$ を求めることで線形回帰をしようとしていたのに、なぜかパラメータが消えてしまいました。代わりに $y$ は平均的には $ f({\bf x}) $ という値になり、ノイズのせいで $ \sigma^2 $ の分散を持ってしまうというちょっと謎のモデルに代わってしまったのです。
そもそも $f({\bf x})$ はどうやって決めるんだよ！？というところがガウシアンプロセスの肝でありますが、ここでは深く突っ込みません。
ガウシアンプロセス回帰では、出力 $y$ は普通の線形回帰と同様ガウス分布に従っているとしており、モデルとして怖気づくほど難しいものではないということを見ました。結局のところ $f({\bf x})$ を求められさえすれば回帰は実行できそうなわけですが、果たして何をどうすれば $f({\bf x})$ を決めることが出来るのでしょうか。
ここで $f({\bf x})$ は（天下り的に非常に不思議なことに）下記のような謎の分布から生成されているということにしてしまいましょう（！？…興味があればPRMLでもSlide Shareでも学べるでしょう）。
$$ f({\bf x}) \sim p(f({\bf x}) \mid {\bf m}, {\bf K}) $$
これは、入力 $\bf x$ に対して、パラメータ $\bf m, K$ を持つ確率分布から、$f({\bf x})$ という値が生起するというようなイメージを持っていただけると良いでしょう。この謎の確率分布のことをガウシアンプロセスと呼び、通常 $\mathcal {GP} ({\bf m}, {\bf K})$ と書いたりします（だいぶいい加減な説明をしてしまっていて、本当はとある $\bf x$ について考えるのではなく、$\bf x_1, ..., x_N$ というように色々な $\bf x_i$ に対して、 $f({\bf x_i})$ という値が生起する同時分布を与えるものが確率過程と呼ばれるものであり、その同時分布がガウス分布であるときにガウス過程という。が細かいことは気にしない）。
兎にも角にもガウシアンプロセス回帰というのは、ガウシアンプロセス $\mathcal {GP} ({\bf m}, {\bf K})$ のパラメータ $\bf m, K$ を推定することが主な仕事になるわけです。
ガウシアンプロセス回帰のまとめ
$$ \begin{align} f({\bf x} ) & \sim \mathcal {GP} ({\bf m}, {\bf K}) \\ y & = f({\bf x}) + \epsilon \\ \epsilon &\sim \mathcal N(0, \sigma ) \end{align} $$
あるいは、キュッとまとめてしまって、あたかも二段階にサンプリングが行われている
$$ \begin{align} f({\bf x} ) & \sim \mathcal {GP} ({\bf m}, {\bf K}) \\ y &\sim \mathcal N(f({\bf x}), \sigma ) \end{align} $$
というモデルであると考えれば良いです。
ガウシアンプロセス回帰の推定
$\bf m, K$（とついでに $\sigma$ ） を推定しようというのが今後の仕事なのですが、こいつらは一体どうやって求めれば良いのでしょうか。 さて、またまた天下りで申し訳ないのですが、カーネル関数 $k(\bf {x_i, x_j})$ というものをパラメトリック関数として選んでしまえば、ガウシアンプロセス回帰の推定は、カーネル関数のパラメータを推定するというものに置き換えることができます。
$\bf m$ も $\bf K $ もカーネル関数を使って表現できてしまうためですが、ここでは具体的な形式は求めません。とにかくカーネル関数なるものを決めてしまえば $\bf m, K$ が決まるということだけ抑えておきましょう。
典型的にはカーネル関数 $k(\bf {x_i, x_j})$ として下記のようなものが選ばれます。
$$ k({\bf x_i, x_j}) = a \exp \left( - \frac{|{\bf x_i - x_j}|^2}{2s^2} \right) $$
こいつは２つのデータ点 $\bf x_i $ と $\bf x_j $ の類似度を表してくれるような関数になっていることが望まれます。 パラメータは $a, s$ です。
まとめ
１．ガウシアンプロセスのモデル
ガウシアンプロセス回帰は $$ \begin{align} f({\bf x} ) & \sim \mathcal {GP} ({\bf m}, {\bf K}) \\ y &\sim \mathcal N(f({\bf x}), \sigma ) \end{align} $$ と書いてやることができる。
２．推定するべきガウシアンプロセスのパラメータ
推定は $\bf m, K$ を求めることである。
３．カーネル関数でガウシアンプロセスのパラメータを書き換える
$\bf m, K$ は ２つのデータ点の類似度を決めるカーネル関数を$k( {\bf x_i, x_j })$ で表現できる
４．推定すべきパラメータをすり替える
カーネル関数を$k( {\bf x_i, x_j })$ をパラメトリックに置く、例えば
$$ k({\bf x_i, x_j}) = a \exp \left( - \frac{|{\bf x_i - x_j}|^2}{2s^2} \right) $$
とすることで、推定は $a, s$ を求めることにすり替わる。
補足
結局のところはカーネル関数なるものがなぜに出てきて、こいつを決めることがガウシアンプロセスを決めることにどう繋がるのか、ということが重要なところなのですが、ここらへんはカーネル法を学ぶことそのものになります。
相応に難しいので詳細は省きますが、カーネル法を駆使するときにはデータ点が重要な役割を担います。なぜなら、カーネル関数は２つのデータ点間の類似度みたいなものなので、２つのデータ点を比べなければカーネル関数の値は求まりませんし、カーネル関数の値が求まらなければガウシアンプロセスも求まらないのです（言い換えれば、カーネル関数のパラメータと更に保持しているデータを使って、ガウシアンプロセスによるサンプリングを計算することになる）。
従って、学習に使った（幾つかの重要な）データ点を保持しなければガウシアンプロセスの計算は行えないということになります。
TensorFlow Probabilityで実践
必要なライブラリのインポート
まずは今回使うライブラリのインポートを行います。
import numpy as np import tensorflow as tf import tensorflow_probability as tfp import matplotlib.pyplot as plt import seaborn as sns sns.set_style() print(tf.__version__) print(tfp.__version__) # -> 1.11.0 # -> 0.4.0
次に頻繁に使うモジュールの名前を下記のように書いておきましょう。
tfd = tfp.distributions psd_kernels = tfp.positive_semidefinite_kernels
これは from tensorflow_probability import distributions as tfd などでも良いのです（個人的にはインポート文が複雑になるのは好きじゃないのでこう書いています）。
でたらめなガウシアンプロセス回帰
データの準備
観測データを作ります。ちなみにこのデータは
$$ y \sim \mathcal N ( (\sin(10x) + \cos(5x)) \exp(-x^2), 0.5) $$
というサンプリングを得ているものです。
observation_noise_variance = .5 def f(x): return (np.sin(10*x[..., 0]) + np.cos(5*x[..., 0])) * np.exp(-x[..., 0]**2) observation_index_points = np.random.uniform(-1., 1., 50)[..., np.newaxis] observations = (f(observation_index_points) + np.random.normal(0., np.sqrt(observation_noise_variance)))
plt.scatter(observation_index_points, observations)
パラメータをフィッティングしていないガウシアンプロセス回帰のサンプリング
ガウシアンプロセス回帰のモデルを適当に作ります。
ガウシアンプロセス（カーネル法）ではカーネル関数の計算にデータ点が必要なため、データ点 $x$ と観測データ $y$ をそれぞれ
observation_index_pointsとobservationsとして与えます。
observation_noise_varianceはガウシアンプロセス回帰の$\sigma^2$に相当する部分です。
index_points = np.linspace(-1., 1., 200)[..., np.newaxis] kernel = psd_kernels.MaternFiveHalves() gprm = tfd.GaussianProcessRegressionModel( kernel=kernel, index_points=index_points, observation_index_points=observation_index_points, observations=observations, observation_noise_variance=observation_noise_variance)
こいつはパラメータが適当に初期化されているので、当然データにフィットしていません。 試しにガウシアンプロセスから $f(x)$ を10個ほどサンプリングしてみましょう。
# ガウシアンプロセスから10回関数fをサンプリングする計算グラフ sample_gprm = gprm.sample(10)
上記の計算グラフを下記で実行します。
with tf.Session() as sess: sample_gprm_ = sess.run(sample_gprm) for i in range(10): plt.plot(index_points, sample_gprm_[i])
もちろんフィッティング前なのでデタラメな関数が１０個表示されているだけです。
学習したガウシアンプロセス回帰
さて、今度はちゃんと学習させたガウシアンプロセス回帰からサンプリングを実行してみましょう。
データ点準備
復習も兼ねて、先程同様の条件でデータ点を作るところからやり直します。
observation_index_points = np.random.uniform(-1., 1., 50)[..., np.newaxis] observations = f(observation_index_points) + np.random.normal(0., .2, 50) plt.scatter(observation_index_points, observations)
ガウシアンプロセスのモデル構築
ガウシアンプロセスの要はパラメトリックなカーネル関数を作ることでした。 今回は下記の数式で表されるカーネル関数を使います。
$$ k({\bf x_i, x_j}) = a \exp \left( - \frac{|{\bf x_i - x_j}|^2}{2s^2} \right) $$
$a$ が amplitudeで、 $2s^2$ がlength_scaleです。いずれも正の値となるようにtf.expで変換しておきます。
amplitude = tf.exp(tf.Variable(np.float64(0)), name='amplitude') length_scale = tf.exp(tf.Variable(np.float64(0)), name='length_scale') kernel = psd_kernels.ExponentiatedQuadratic(amplitude, length_scale)
こうすることで、
$$ \begin{align} f({\bf x} ) & \sim \mathcal {GP} ({\bf m}, {\bf K}) \\ y &\sim \mathcal N(f({\bf x}), \sigma ) \end{align} $$
における１つ目の式のパラメータを定式化できたことになります（$\bf m$ も $\bf K$ もカーネル関数を決めれば表現できるのでした）。 ついでに２つ目の式の $\sigma$ も下記のコードで表しておきましょう。
observation_noise_variance = tf.exp(
tf.Variable(np.float64(-5)), name='observation_noise_variance')
こいつらをガウシアンプロセスを表すtfd.GaussianProcessに食わせてやれば、モデル全体を定式化できたことになります。
gp = tfd.GaussianProcess(
kernel=kernel,
index_points=observation_index_points,
observation_noise_variance=observation_noise_variance)
損失関数の設定
最尤推定をするためには、仮定したモデルにおけるデータ生成の対数尤度関数を作り、それを最大化すればいいです。
TensorFlowでは典型的に最小化問題として表すのが通例であるので、対数尤度関数の符号を変えたものを損失関数とみなして最小化を行うことで最尤推定を実行します。
neg_log_likelihood = -gp.log_prob(observations) optimizer = tf.train.AdamOptimizer(learning_rate=.05, beta1=.5, beta2=.99) optimize = optimizer.minimize(neg_log_likelihood)
対数尤度を求めるのはTensorFlow Probabilityでは極めて簡単で、確率分布のクラスがlog_probメソッドを持っており、観測データを与えてやるだけで対数尤度を返してくれます。コードではgp.log_prob(observations)がそれに該当します。
あとは慣れたようにoptimizerを作って、minimizeメソッドに損失関数を渡すことで最小化のための計算グラフが構築できます。
ガウシアンプロセス回帰のサンプラー
gp = tfd.GaussianProcess()のインスタンスにkernelを持たせて、最尤推定する計算グラフは作れたので、後は同じkernelを参照してサンプリングを行うtfd.GaussianProcessRegressionModelの計算グラフを作っておきます。
index_points = np.linspace(-1.5, 1.5, 200)[..., np.newaxis] gprm = tfd.GaussianProcessRegressionModel( kernel=kernel, index_points=index_points, observation_index_points=observation_index_points, observations=observations, observation_noise_variance=observation_noise_variance) samples = gprm.sample(10)
最初の実験同様、関数のサンプリングを１０回行う計算グラフを作っておきましょう。 （ココらへんの書き方は完全にDefine and RunであるTensorFlowによるものなので、慣れが必要です）
いざ学習！
準備万端なので、実際に学習を行います。
with tf.Session() as sess: sess.run(tf.global_variables_initializer()) for i in range(1000): _, neg_log_likelihood_ = sess.run([optimize, neg_log_likelihood]) if i % 100 == 0: print("Step {}: NLL = {}".format(i, neg_log_likelihood_)) print("Final NLL = {}".format(neg_log_likelihood_)) samples_ = sess.run(samples) plt.scatter(np.squeeze(observation_index_points), observations) plt.plot(np.stack([index_points[:, 0]]*10).T, samples_.T, c='r', alpha=.2)
Final NLL = 10.56364360237319となり、ガウシアンプロセス回帰による１０個の関数のサンプリングは下記となりました。
更に進むために
ガウシアンプロセスをもっとベイズ的に取り扱う場合にはカーネル関数のパラメータに事前分布を置くという方法があります。この場合はMCMCを活用することでTensorFlow Probabilityで簡単に計算が実行できるのですが、イマイチまだ使い方がよく分かってないので、パスで…。
