読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

情報量規準LOOCVとWAICの比較

この記事はStan Advent Calendar 2016およびR Advent Calendar 2016の12月7日の記事です。StanコードとRコードは記事の最後にあります。

背景は以下です。

  • [1] Aki Vehtari, Andrew Gelman, Jonah Gabry (2015). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. arXiv:1507.04544. (url)
  • [2] 渡辺澄夫. 広く使える情報量規準(WAIC)の続き (注4)【WAICとクロスバリデーションの違いについて】 (url)
  • [3] Sumio Watanabe. Comparison of PSIS Cross Validation with WAIC. (url)

leave-one-outクロスバリデーション(以下LOOCV)およびWAICは予測のよさをベースにしたモデル選択に用いられる情報量規準であり、ともに汎化誤差(Generalization Error、以下GE)の近似です。それにもかかわらず[1]では、本来性能評価では必須と思われる汎化誤差との比較がありません。実データ(真のモデルが未知の状況)で用いているためと思いますが、これではいけないように思います。この記事では僕が日常的に使用するような5つの基本的なモデルを使い、真のモデルが既知の状況でGE・LOOCV・WAICの性能比較を行いました。

具体的には[2]に以下のようなコメントがあります。なお[1]ではPareto Smoothed Importance SamplingでLOOCVを算出しており、PSISCVとも呼ばれるようです。

  • (0) まず同じデータに対してマルコフ連鎖モンテカルロ法を何度も行ったときの値の揺らぎを調べてみましょう.WAICの分散はISCVおよびPSISCVの分散よりも小さくなります。つまり、WAICはISCVおよびPSISCVよりもマルコフ連鎖揺らぎに対して強いということができます.
  • (2) しかしながら,CVもWAICも汎化誤差を推定することが本来の目的です.CVとWAICの両方の厳密値が計算できたとして,(つまりMCMC法で無限にサンプルが取れたとき),どちらの方が汎化誤差の推定として優れているのでしょうか。(中略) 我々の実験では,GEを汎化誤差とするとき,ほぼ,いつでも E[|PSISCV-GE|] > E[|WAIC-GE|] が成り立つのですが・・・。このページをご覧の皆様にはぜひ実験してみていただければと思います。なお、E[ ]は学習用データのでかたについての平均を表しています。

そこでこれら2点について検証してみました。先に「(2)のコメント」について説明します。

GE・LOOCV・WAICの比較

使用した5つのモデルとシミュレーションの手順を説明します。

重回帰

真のモデルは以下です。

  Y \sim Normal(1.3 - 3.1 X_1 + 0.7 X_2, 2.5)

あてはめたモデルは以下です。

  Y \sim Normal(b_1 + b_2 X_1 + b_3 X_2, \sigma)

データ点の数Nについては10,30,100,300を試しました。例としてN = 10の場合を説明します。まず乱数でデータX(すなわち X_1, X_2)を生成します。次にそのXの値を使って1000通りのYを生成します(学習用データのでかたの平均をとるため)。各YについてStanでiter=11000, warmup=1000, chains=4で実行して合計40000個のMCMCサンプルを得てGE・LOOCV・WAICと「LOOCV - GE」と「WAIC - GE」を求めました。その後、Nごとに「LOOCV - GE」と「WAIC - GE」のboxplotを描きました。

ロジスティック回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

  Y \sim Bernoulli(inv\_logit(0.8 - 1.1 X_1 + 0.1 X_2))

あてはめたモデルは以下です。

  Y \sim Bernoulli(inv\_logit(b_1 + b_2 X_1 + b_3 X_2))

  b_1,b_2,b_3 \sim Student\_t(4,0,1)

非線形回帰 ミカエリス・メンテン型

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

  Y \sim Normal(10.0 X / (2.0 + X), 0.8)

あてはめたモデルは以下です。

  Y \sim Normal(m X / (k + X), \sigma)

  k \sim Uniform(0, 12)

  m \sim Uniform(0, 20)

真のモデルが含まれない場合

あてはめたモデルが以下の場合も試しました。

  Y \sim Normal(a + b X, \sigma)

ノイズがt分布に従う重回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

  Y \sim Student\_t(4, 1.3 - 3.1 X_1 + 0.7 X_2, 2.5)

あてはめたモデルは以下です。

  Y \sim Student\_t(4, b_1 + b_2 X_1 + b_3 X_2, \sigma)

階層モデル

手順は重回帰の場合とおおよそ同じですが、データ点の数とモデルが異なります。グループの数を10に固定し、データ点の数Nについては20,50,130,400を試しました(それぞれ各グループで2,5,13,40人)。また真のモデルは以下です。

  \mu\left[g\right] \sim Normal(0, 10)   g = 1,2,\dots,10

  Y\left[n\right] \sim Normal(\mu\left[g\left[n\right]\right], 2.0)   n = 1,2,\dots,N

あてはめたモデルは以下です。

  \mu\left[g\right] \sim Normal(\mu_0, \sigma_0)   g = 1,2,\dots,10

  Y\left[n\right] \sim Normal(\mu\left[g\left[n\right]\right], \sigma)   n = 1,2,\dots,N

  \mu_0,\sigma_0 \sim Student\_t(4, 0, 10)

階層モデルにおいては何を予測したいのか(どういう状況の汎化誤差を知りたいのか)を注意深く考える必要があります。以下の記事を参照。

statmodeling.hatenablog.com

ここでは以下の2つの場合を計算しました。

  • 既存の各グループに、新しく1人ずつ加わる場合
  • 別の新しいクラスができて、新しく1人が加わる場合

結果

重回帰

大きな差はありませんでした。Nが小さい場合にWAICの方がわずかにGEに近くなる傾向があるようです。

ロジスティック回帰

大きな差はありませんでした。Nが小さい場合にLOOCVの方がわずかにGEに近くなる傾向があるようです。

非線形回帰 ミカエリス・メンテン型

大きな差はありませんでした。Nが小さい場合にWAICの方がわずかにGEに近くなる傾向があるようです。

  • 真のモデルが含まれない場合

この場合はNを増やしても汎化誤差引くエントロピー(=予測分布と真の分布のKL情報量)は0に近づかず、0.89ほどで下げ止まります。そして、Nが小さい場合にWAICの方がGEに近くなる傾向があるようです。

ノイズがt分布に従う重回帰

Nが小さい場合にWAICの方がGEに近くなる傾向があるようです。なおboxからはみ出るoutlierの値が大きく、そのままプロットすると見づらくなるので、図の縦軸を制限しました。

階層モデル

  • 既存の各グループに、新しく1人ずつ加わる場合

Nが小さい場合にWAICの方がGEに近くなる傾向があるようです。

  • 別の新しいクラスができて、新しく1人が加わる場合

グループ数が少ない場合、グループあたりの人数を増やしてもLOOCV - GEおよびWAIC - GEの中央値は0に近づいていきません。グループ数が少ないと、グループあたりの人数を増やしてもグループを生成する平均と標準偏差のパラメータは精度よく求められないことを反映しているのだと思います。全体的にWAICの方がわずかにGEに近くなる傾向があるようです。

階層モデル その2 2016.12.14追記

伊庭先生から以下のような要望がありました。

そこでモデルは同じものを用いて、グループあたりの人数を固定し(2,5,13のうちいずれか)、グループの数が10,30,100の各場合で実行してみました。

  • 既存の各グループに、新しく1人ずつ加わる場合

グループあたりの人数が少ない場合、グループ数を増やすとLOOCV - GEおよびWAIC - GEのバラツキは少なくなるもの、それらの中央値は0に近づいていきません。グループ数を増やしても各グループには2人しかいないため、グループごとの予測はあたらないままということを反映しているのだと思います。グループあたりの人数が小さい場合にWAICの方がGEに近くなる傾向があるようです。

  • 別の新しいクラスができて、新しく1人が加わる場合

Gが小さい場合にWAICの方がわずかにGEに近くなる傾向があるようです。

* * *

5つのモデルを通して見ると、Nが小さい場合、すなわち1サンプルの重みが大きい場合にはWAICのほうがLOOCVよりもよい汎化誤差の近似になっているようです。[2]の「PSISCVとWAIC:実験例 追加」によると、影響の大きい(重みの大きい)サンプルが存在する場合にも同じような結果になるようです。以下の伊庭先生のツイートはこのような状況を指していると思われます。

WAICはLOOCVよりMCMCの揺らぎに強いか?

次に「(0)のコメント」について検証しました。前述の5つのモデルに対し、N10または100とし、各Nについてデータを5通り生成しました。さらに各データセットに対して、MCMCのシードを1000通り試し、1000個のLOOCVとWAICを求め、それぞれの分散と平均と変動係数を求めました。

結果

いずれの場合についても実用上の差が認められませんでした。N100の場合より10の場合の方がMCMCの揺らぎがありますが、それでも変動係数にして高々1%程度でした。なお、この結果はMCMCサンプルを求めるアルゴリズムの違いやMCMCサンプルの数にも依存すると思います。

まとめ

  • Nが大きい場合には、WAICとLOOCVにほぼ差がないと言えるでしょう。
  • Nが小さい場合には、WAICの方が汎化誤差のよい近似になっていると言えるでしょう。
  • さらに理論的な美しさや計算速度も含めて総合的に判断すると、WAICに軍配が上がると思います。

おまけ

データ Xとは異なる xにおける yの予測分布 y_{pred}(y|x)を考えたい場合があると思います。その場合は、一般によい情報量規準があるか未解明で、研究対象として興味深いようです(渡辺澄夫先生(私信))。例えば単回帰の場合には、真のモデルがあてはめたモデルに含まれており、かつモデルに依存する量を使うと情報量規準に準ずる量を構築できるようです。興味深いです。


ソースコード

「(2)のコメント」を検証するための、重回帰の場合と階層モデルの場合のStanコードとRコードを公開します。

重回帰

Stanコード

重回帰のモデルの部分は「StanとRでベイズ統計モデリング」の9.2.1項と同じです。異なる点はgenerated quantitiesブロックでGE・LOOCV・WAICを算出するためにlog_liky_predを算出している点です。

Rコード

例としてN100の場合を載せます。理解しやすさのため、並列化していないコードにしてありますが、実際には色々並列化して計算しています。

  • 16~17行目: 学習用データの出かたの平均をとりたいので、データYを乱数で生成しています。
  • 22~31行目: データごとのエントロピーと汎化誤差を計算しています。以前の記事参照。このコードのように十分なMCMCサンプルから予測分布の近似関数を求めて汎化誤差を算出する方法のほか、直接MCMCサンプルを使って求める方法もあるかと思います(渡辺先生はそうしています)。
  • 23~24行目: 予測分布は滑らかだと仮定し、40000個のMCMCサンプルから予測分布の密度関数を計算しています。Rのデフォルトのdensity関数よりも{KernSmooth}パッケージのbkde関数の方が優秀っぽいのでこちらを使っています(参考pdf)。
  • 25行目: 真の分布です。
  • 26行目: エントロピーの計算の際に積分される関数です。
  • 27行目: 汎化誤差の計算の際に積分される関数です。予測分布の密度推定した結果をapproxfunで関数に変換している関係上、f_predがごくまれに絶対値の小さな0以下の値を返します。それを避けるためにifelse関数をかませてあります。
  • 28~29行目: それぞれエントロピーと汎化誤差を計算しています。ここではf_true正規分布なので、-6SD~+6SDまで積分すれば十分よい近似となります。Rのデフォルトのintegrate関数はちょっと賢くてadaptiveに積分しているので計算は早いのですが、まれに不安定で計算ができない場合があります。そのため、少し遅いですが安定な数値積分の手法であるRombergの方法を使っています({pracma}パッケージのromberg関数または{Rmpfr}パッケージのintegrateR関数を使うことができます)。
  • 34行目: 汎化誤差(GE)のサンプルに関する平均を求めています。学習用データと同じNXの値を持つ新しいデータセットに対して予測を行い、1サンプルあたりの汎化誤差を求めていることに相当します。
  • 35~36行目: Stanのチームが開発している{loo}パッケージを用いてlooicwaicを求めています。彼らの情報量規準のスケールはAICやDICとあわせて 2N倍となっているので、1サンプルあたりにスケールをあわせる意味で2*Nで割っています。

階層モデル「既存の各グループに、新しく1人ずつ加わる場合」

階層モデルの場合のWAICの詳しい解説は以前の記事を参照してください。

Stanコード

Rコード

重回帰の場合と似ています。

  • 23~32行目: グループgに1人加わった場合の汎化誤差error_by_groupを求めています。重回帰の場合にサンプルごと(nごと)だったのが、グループごと(gごと)にインデックスが変わっただけで、内容は変わっていません。
  • 35行目: 「既存の各グループに、新しく1人ずつ加わる場合」なので、各グループの汎化誤差の和になります。
  • 36~37行目: 各グループごとにLOOCVまたはWAICを求めて和をとっています。

階層モデル「別の新しいクラスができて、新しく1人が加わる場合」

Stanコード

Rコード

  • 25~26行目: この場合は真の分布が積分を含んでいるので少し複雑です。

残りは「既存の各グループに、新しく1人ずつ加わる場合」とほぼ同じです。

なお、この記事は以下のツイートによってやってみようかなと思いました。

GPy(Pythonのガウス過程用ライブラリ)の使い方

概要

GPyを用いて、サンプルパスの生成、ガウス過程回帰、クラス分類、ポアソン回帰、Bayesian GPLVMを実行しました。自分用のメモです。

参考資料

理論的背景は上記の[3]を参考にしてください。日本語でもガウス過程の解説がMLPシリーズから豪華著者陣で出るようです。超期待しています。

以下のサンプルプログラムは基本的に[2]を元にしています。しかし、古くてそのままでは動かないプログラムや分かりにくいプログラムを少し加工修正しています。なお、環境は以下の通りです。

サンプルパスの生成

RBFカーネルで適当に定めたパラメータの値でサンプルパスを生成するプログラムです。カーネルそのものやカーネルのパラメータを変えることでどのようなサンプルパスを生成するのかシミュレーションしたい場合によく使います。

import GPy
import numpy as np
import matplotlib.pyplot as plt

kernel = GPy.kern.RBF(input_dim=1, variance=1, lengthscale=0.2)

np.random.seed(seed=123)
N_sim = 100
x_sim = np.linspace(-1, 1, N_sim)
x_sim = x_sim[:, None]
mu = np.zeros(N_sim)
cov = kernel.K(x_sim, x_sim)
y_sim = np.random.multivariate_normal(mu, cov, size=20)

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
for i in range(20):
    ax.plot(x_sim[:], y_sim[i,:])
fig.savefig('output/fig1.png')
  • 5行目: これでカーネルを定義します。入力の次元(input_dim)は必須です。
  • 10行目: GPyの関数の多くは、引数のshape(データ点の数, 1)である必要があります。そこで[:, None]を加えてその形にしています。
  • 12行目: kernelオブジェクトに対しK関数を使うと分散共分散行列を作成できます。
  • 13行目: numpyの関数で多変量正規分布からサンプルを生成しています。

ガウス過程回帰(入力1次元・出力1次元)

手順は カーネルを定める→モデル作成→最適化 だけです。

import GPy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

kernel = GPy.kern.RBF(1)
# kernel = GPy.kern.RBF(1) + GPy.kern.Bias(1) + GPy.kern.Linear(1)

d = pd.read_csv('http://kasugano.sakura.ne.jp/images/2016/20161112/data-GPbook-Fig2_05.txt')
model = GPy.models.GPRegression(d.X[:, None], d.Y[:, None], kernel=kernel)
model.optimize()
model.plot()
plt.savefig('output/fig2.png')

## prediction
x_pred = np.linspace(-10, 10, 100)
x_pred = x_pred[:, None]
y_qua_pred = model.predict_quantiles(x_pred, quantiles=(2.5, 50, 97.5))[0]
  • 6行目: RBFカーネルinput_dim = 1で作成しています。7行目はRBF+Bias+Linearのカーネルを使う場合です。足したり掛けたりするだけで複雑なカーネルを作ることができるインターフェースが素敵です。
  • 10行目: GPy.models.GPRegression関数でモデルを作成しています。print(model)m['']とするとモデルに含まれるパラメータを見ることができます。特に指定しなければ、すべてのパラメータが最適化の対象となります。ちなみに各パラメータに固定値を与えることや制限をかけることができます(詳しくはこれこれを参照)。
  • 11行目: 最適化をしています。オプションでiterationの数などを指定できます。
  • 12~13行目: 最適化後のモデルをプロットしています。
  • 16~18行目: 最適化後のモデルを使って予測を行っています。

ガウス過程回帰(入力2次元・出力1次元)

import GPy
import numpy as np
import matplotlib.pyplot as plt

kernel = GPy.kern.Matern52(2, ARD=True)

np.random.seed(seed=123)
N = 50
X = np.random.uniform(-3.,3.,(N, 2))
Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2]) + np.random.randn(N,1)*0.05

model = GPy.models.GPRegression(X, Y, kernel)
model.optimize(messages=True, max_iters=1e5)
model.plot()
plt.savefig('output/fig3.png')

model.plot(fixed_inputs=[(0, -1.0)], plot_data=False)
plt.savefig('output/fig3-slice.png')

## prediction
x_pred = np.array([np.linspace(-3, 3, 100), np.linspace(3, -3, 100)]).T
y_qua_pred = model.predict_quantiles(x_pred, quantiles=(2.5, 50, 97.5))[0]
  • 5行目: 今回はMatern5/2カーネルを使っています。オプションのARD=Trueは入力の次元1つに対し、1つのlengthscaleパラメータを割り振ること(すなわちGPは等方でないことを表します)。
  • 17~18行目: 2次元の入力のうち一部を固定した図(スライスした図;2枚目の図)を描いています。ここでは、fixed_inputs=[(0, -1.0)]でインデックス0の入力を-1.0に固定しています。
  • 21~22行目: 入力1次元のガウス過程回帰と同様に予測をしています。入力の次元に注意です。

スパースなガウス過程回帰

補助変数法やコンパクトなガウス過程回帰とも呼ばれます。ガウス過程はデータ点の数N逆行列を求める必要があり、その部分にN^3のオーダーの時間がかかります。そのため、データ点が増えると次第に遅くなります。そこで、一部の補助変数(inducing inputs)を入力次元の代表点として扱い、対数尤度を近似することで計算を高速化させる方法があります。それがこの節の方法になります。

import GPy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

kernel = GPy.kern.RBF(1)

d = pd.read_csv('http://kasugano.sakura.ne.jp/images/2016/20161112/data-GPbook-Fig2_05.txt')
m_full = GPy.models.GPRegression(d.X[:, None], d.Y[:, None], kernel=kernel)
m_full.optimize()

Z = np.hstack((np.linspace(-6, -3, 3), np.linspace(3, 6, 3)))[:,None]
# Z = np.linspace(-6, 6, 12)[:, None]
m_sparse = GPy.models.SparseGPRegression(d.X[:, None], d.Y[:, None], Z=Z)
m_sparse.optimize()
m_sparse.plot()
plt.savefig('output/fig4.png')
print(m_sparse.log_likelihood(), m_full.log_likelihood())
  • 12行目: 補助変数の初期値です。6個をテキトーに定めました。結果は1枚目の図です。
  • 13行目: こちらは12個の場合です。結果は2枚目の図です。
  • 14行目: GPy.models.SparseGPRegression関数を使います。補助変数はZで指定します。
  • 15行目: 補助変数の位置も最適化の対象となります。
  • 18行目: modelオブジェクトに対しlog_likelihood関数を使うと対数尤度を取得できます。最適化の後の対数尤度を見ると、補助変数6個の場合が-28.85、補助変数12個の場合が-18.02、補助変数を使わないフルモデルの場合が-17.92となりました。12個の補助変数で十分よく近似できていることが分かります。

クラス分類

ここではPRML下のFig.6.12相当の図を再現してみます。

import GPy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

kernel = GPy.kern.RBF(2, ARD=True)

d = pd.read_csv('http://kasugano.sakura.ne.jp/images/2016/20161112/data-classification.txt')
model = GPy.models.GPClassification(d[['X1', 'X2']].values, d.Y[:, None])
model.optimize()

ax = model.plot(plot_data=False)
d0 = d[d.Y == 0]
d1 = d[d.Y == 1]
ax.plot(d0.X1, d0.X2, 'ro')
ax.plot(d1.X1, d1.X2, 'bo')
plt.savefig('output/fig5.png')
  • 9行目: GPy.models.GPClassification関数でクラス分類のモデルを組み立てることができます。

ポアソン回帰

久保緑本の11章の欠測値なしのモデルを実行します。

import GPy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

kernel = GPy.kern.RBF(1)

d = pd.read_csv('http://kasugano.sakura.ne.jp/images/2016/20161112/data-kubo11a.txt')

model = GPy.core.GP(X=np.linspace(1, 50, 50)[:,None], Y=d.Y[:,None], kernel=kernel,
    inference_method=GPy.inference.latent_function_inference.Laplace(),
    likelihood=GPy.likelihoods.Poisson())
model.optimize()
model.plot()
plt.savefig('output/fig6.png')

x_pred = np.linspace(1, 50, 50)[:, None]
f_mean, f_var = model._raw_predict(x_pred) # Predictive GP for log intensity mean and variance
f_upper, f_lower = f_mean + 2*np.sqrt(f_var), f_mean - 2.*np.sqrt(f_var)
plt.plot(x_pred, np.exp(f_mean), color='blue', lw=2)
plt.fill_between(x_pred[:,0], np.exp(f_lower[:,0]), np.exp(f_upper[:,0]), color='blue', alpha=.1)
plt.savefig('output/fig6-mean.png')
  • 10~12行目: 少し凝ったモデルを使用したい場合には、GPy.core.GP関数を使って尤度を自分で設定する必要があります。推定方法もあわせて指定します。ここでは単純なポアソン回帰なので、用意されているGPy.likelihoods.Poisson関数を使えば完了です。
  • 17~21行目: 真の平均の推定値と±2SDのグラフ(2枚目の図)を描いています。

ガウス過程回帰(入力2次元・出力2次元)

出力が2次元となると、モデルの選択肢が増えます。その前に「どうして複数次元の出力が必要なのか?各出力は相関しているのか(一方の出力が他方の出力を予測するヒントになるのか)?」といった問いが重要だとNeil Lawrenceは述べています(メーリングリストより)。もしそれらの問いの回答がNoならば、出力1次元のモデルを複数組み合わせたモデルの方がよいかもしれません。

ここでは、「出力間に相関がある簡単なモデル」「出力間に相関がないモデル」「出力間に相関がある凝ったモデル」の順にすすめます。

import GPy
import numpy as np
import matplotlib.pyplot as plt

f_output1 = lambda x: 4*np.cos(x/5) - 0.4*x - 35 + np.random.rand(x.size) * 2
f_output2 = lambda x: 6*np.cos(x/5) + 0.2*x + 35 + np.random.rand(x.size) * 8

np.random.seed(seed=123)
X1 = np.random.rand(100)
X2 = np.random.rand(100)
X1 = X1*75
X2 = X2*70 + 30
Y1 = f_output1(X1)
Y2 = f_output2(X2)

x_pred1 = np.random.rand(100)*100
x_pred2 = np.random.rand(100)*100
y_pred1 = f_output1(x_pred1)
y_pred2 = f_output2(x_pred2)

def plot_2outputs(m):
    fig = plt.figure(figsize=(12, 8))
    ax1 = fig.add_subplot(211)
    ax1.set_ylim([-120, -20])
    ax1.set_title('Output 1')
    m.plot(plot_limits=[0, 100], fixed_inputs=[(1, 0)], which_data_rows=slice(0, 100), ax=ax1)
    ax1.plot(x_pred1, y_pred1, 'rx', mew=1.5)
    ax2 = fig.add_subplot(212)
    ax2.set_ylim([-20, 100])
    ax2.set_title('Output 2')
    m.plot(plot_limits=[0, 100], fixed_inputs=[(1, 1)], which_data_rows=slice(100, 200), ax=ax2)
    ax2.plot(x_pred2, y_pred2, 'rx', mew=1.5)


K = GPy.kern.Matern32(1)
B = GPy.kern.Coregionalize(input_dim=1, output_dim=2)
kernel = K**B
model = GPy.models.GPCoregionalizedRegression(X_list=[X1[:, None],X2[:, None]], Y_list=[Y1[:, None],Y2[:, None]], kernel=kernel)

model['.*Mat32.var'].constrain_fixed(1)
model.optimize()
plot_2outputs(model)
plt.savefig('output/fig7a.png')

model['.*coregion.W'].constrain_fixed(0)
model.randomize()
model.optimize()
plot_2outputs(model)
plt.savefig('output/fig7b.png')


K1 = GPy.kern.RBF(1)
K2 = GPy.kern.Bias(1) + GPy.kern.Linear(1)
B1 = GPy.kern.Coregionalize(1, output_dim=2)
B2 = GPy.kern.Coregionalize(1, output_dim=2)
kernel = K1**B1 + K2**B2

X = np.vstack((np.concatenate([X1, X2]), np.hstack((np.zeros(100), np.ones(100))))).T
Y = np.hstack((Y1, Y2))[:, None]
model = GPy.models.GPRegression(X, Y, kernel)
model.optimize()
plot_2outputs(model)
plt.savefig('output/fig7c.png')


## prediction
x_pred = np.arange(100, 110)[:, None]
x_pred = np.hstack([x_pred, np.ones_like(x_pred)])
output_index_pred = {'output_index':x_pred[:,1:].astype(int)}
y_pred = model.predict(x_pred, Y_metadata=output_index_pred)
  • 5~19行目: デモデータ作成部分です。
  • 21~32行目: プロットする関数を定義しています。
  • 36行目: 出力間の関係を定める行列B(coregionalization matrix)を作成しています。詳しくはここを参照。
  • 37行目: kernelオブジェクトに対する**クロネッカー積となります。
  • 38行目: GPy.models.GPCoregionalizedRegression関数を使うことで、複数の出力の次元が相関を持ち、出力の各次元でノイズの大きさが異なるモデルを簡単に作成することができます。なお、入力と出力はlistで渡します。
  • 40行目: 38行目でもモデルにノイズが入るので、Matern3/2カーネルに含まれるノイズを1に固定しています。
  • 41~43行目: このモデルの結果は1枚目の図です。
  • 45行目: Bを対角行列に固定しています。出力次元ごとにガウス過程回帰を行うのと同じになります。
  • 46行目: 41行目で最適化された値になっているので、いったん初期値をぐちゃぐちゃにする意味です。
  • 47~49行目: このモデルの結果は2枚目の図です。1枚目の図との違いに注目してください。
  • 52~56行目: coregionalization matrixをカーネルの種類ごとに用意して組み合わせることもできます。
  • 60行目: GPy.models.GPRegression関数を使うと、すべての出力次元の共通の大きさのノイズとなります。なお、GPy.models.GPRegression関数を用いる場合にはXは「1列目に値・2列目に出力次元のインデックス」となっているndarrayを渡す必要があります。YXに対応するndarrayです。
  • 61~63行目: このモデルの結果は3枚目の図です。予測区間が狭いです。モデルが複雑で過学習の恐れがあるかもしれません。
  • 67~70行目: 予測の例です。出力が複数ある場合にはdictionaryを作って渡します。ここではoutput_index1(すなわち2番目の出力)が100110でどのような出力になるか予測しています。

Bayesian GPLVM

前の記事と同じようにPRMLでおなじみのOil Flowのデータに対してBayesian GPLVMを実行します。

from scipy.io import loadmat
import scipy.io as spio
import GPy
import matplotlib.pyplot as plt

d = spio.loadmat('input/3Class.mat')
X = d['DataTrn']
X -= X.mean(0)
L = d['DataTrnLbls'].nonzero()[1]
input_dim = 2 # How many latent dimensions to use

kernel = GPy.kern.RBF(input_dim, ARD=True) + GPy.kern.Bias(input_dim) + GPy.kern.Linear(input_dim) + GPy.kern.White(input_dim)
model = GPy.models.BayesianGPLVM(X, input_dim, kernel=kernel, num_inducing=30)
model.optimize(messages=True, max_iters=5e3)
model.plot_latent(labels=L)
plt.savefig('output/fig8.png')
  • 10行目: 潜在変数の次元。ここではチュートリアルと同じように2にしています。
  • 13行目: GPy.models.BayesianGPLVM関数の一発で補助変数込みでモデルが作成できます。

不明点

GPy.models.GPRegression関数を使うと、パラメータにノイズの大きさが加わります。このノイズと、カーネル側でGPy.kern.White関数で設定したノイズの違いがよく分かりません。なお、簡単なモデルで両方とも設定すると数式の上では識別できなくなると思うのですが、最適化の結果は分散をちょうど半分に分ける形で推定されて全体的な予測分布は変わりません。

Bayesian GPLVMをStanで実装してみた

この記事の続きです。PRML下の12章に出てくるOil Flowのデータ(データ点1000個×特徴量12個)に対してBayesian GPLVMで2次元(または3次元)の潜在変数空間にマッピングして綺麗に分離されるか見てみます。

まずはPRMLにもあるように普通の主成分分析でやると以下になります。綺麗には分離されません。

次にBayesian GPLVMでやってみます。Stanコードは以下になります。

  • 2~4行目: NKDはそれぞれ、データ点の数・特徴量の数・最終的に落とし込む潜在空間の次元です。
  • 14行目: 潜在変数です。
  • 15行目: カーネルに含まれるパラメータです。僕が理解したところだと特徴量ごとにガウス過程が存在するのでKごとに異なる値を持つようにしています。
  • 19行目: 同様に特徴量ごとに分散共分散行列があります。
  • 20・22~23行目: カーネルの定義で効率的な行列演算をするため、matrix型をvector型の配列に持ち替えます。
  • 24~28行目: カーネルの定義です。ここではGaussian+bias+linear+white noiseにしました。カーネルについてはGP summser school 2015のKernel Designの講義資料 (pdf)The Kernel Cookbookなどを参照してください。
  • 29行目: 潜在空間に対する縛りです。代わりにparametersブロックでlowerupperを定めてもOKです。
  • 30~31行目: カーネルに含まれるパラメータの事前分布です。軽くしばっています。ある範囲の一様分布にするなど他にも色々考えられると思います。
  • 32~33行目: ガウス過程の部分です。この書き方だと各特徴量は独立になっていますが、さらに特徴量間の相関を考慮したモデル(例えばCoregionalized Regression Model)もあります。ここではすでに計算量が膨大なので独立としました。

Rからの実行方法は以下です。計算が重いのでADVIを使いました。

  • 5行目: データがmatlabのファイルで与えられていたのでそれを読み込んでいます。結果はリストになります。
  • 8行目: 潜在変数の次元Dを与えています。最終的に2次元にプロットしたいです。余裕をもって次元を設定し、主成分分析のように寄与の大きい次元トップ2だけを抽出することで2次元に射影する方法もあるようです。ここでははじめから2次元の空間とします。また3次元の空間も試して3Dプロットしてみます。
  • 10~16行目: まずはふつうの主成分分析しています。前述の図はここで出力しているresult-pca.pngになります。またGPLVMを実行する際の初期値にする意味もあります。
  • 18~28行目: Bayesian GPLVMを実行しています。21行目ではinitオプションでPCAの結果を初期値として設定しています。22行目では時間短縮のためetaを求めないで1で与えています(手元のモデルとデータではチューニングの結果がいつもeta = 1だったという理由があります)。

計算時間はAWS EC2のc4.xlargeを使っても3時間半ぐらいかかりました。かなり遅いです。ADVIの代わりにStanのoptimizing関数を使って推定した方がよいかもしれません。また汎用的な確率的プログラミング言語ではガウス過程に特化した専用ライブラリにはかないません。Stanのモデルはユーザの問題にあわせた拡張が簡単なので、その点で使う価値はあると思います。特にモデルを拡張する予定がないならば、もしくはデータが巨大ならば、Pythonガウス過程に特化したライブラリであるGPyなどの使い方を学ぶべきと思います。

結果

D = 2の場合

潜在変数xの乱数サンプルの中央値を使って2次元の散布図を描くと以下になりました(result-bgplvm.png)。それなりに綺麗に分離していると思います。乱数の種の影響も見ましたが、おおよそ似たような結果になりました。

D = 3の場合

同様に3次元の散布図を描きました。下の図はマウスでグリグリできます。2次元の場合より若干綺麗に分離していそうです。