情報論的学習理論テクニカルレポート 2009
Technical Report on Information-Based Induction Sciences 2009 (IBIS2009)

カーネルマルコフ連鎖モンテカルロ法による測定誤差モデル推定

Measurement Error Model Estimation byKernel Markov Chain Monte Carlo Method
赤穂 昭太郎∗Shotaro Akaho

伊庭 幸人†
Yukito Iba

Abstract: Measurement error models are the statistical models in which random noise isadded to input variables as well as output variables. We consider the estimation problem of
the regression function in a reproducing kernel Hilbert space (RKHS) for the measurement
error models. We apply Markov chain Monte Carlo approach to estimate the posterior ofthe function. To deal with the infinite dimensionality of RKHS, we introduce a trick to
exchange the order of sampling of the hidden variable and the function.
Keywords: Bayesian framework, Markov chain Monte Carlo, reproducing kernel Hilbert
space, Gaussian process, measurement error model

1

まえがき

基底ベクトルをもつカーネル関数を用いれば無限次元の入力値から出力値への関数を推定する回帰問題では，通常入力値の関数にノイズが加わったものが観測されていると仮定することが多い．ところが，位相応答曲線の推定 [1, 2] など，一般には入力値そのものにも観測ノイズが加わっている場合も多い． すなわち，

x = z + x ,y = f (z) + y

自由度をもつ関数空間となるので，あらかじめ関数クラスが規定できないような場合には柔軟な関数あてはめができる可能性がある．
無限自由度の関数によるあてはめを実現するためにはなんらかのトリックが必要である．出力にしかノイズが乗っていない通常の回帰で最適解 (MAP 推定値) のみが(1)必要な場合は、観測点におけるカーネル関数だけを使って解析的に解が求められる (Representer 定理)．マルコのように潜在変数 z を介して入力 x にも出力 y にもフ連鎖モンテカルロ法 (MCMC) などを用いて事後分布ノイズがのる．このようなモデルは「測定誤差モデル」からのサンプルを生成する場合には、サンプルそのもの(Measurement Error Model) として知られている [3]．

は無限自由度であって直接に扱えない．しかし，関数値関数 f のクラスとしては直線 [4] やスプライン関数を評価する点の集合があらかじめ与えられた場合には，[5] の場合が研究されているが，本稿では再生核ヒルベそれに対応する基底を加えることで，有限次元の範囲でルト空間 (RKHS: Reproducing Kernel Hilbert Space)処理することができる [8, 9]．

上の関数を考える．再生核ヒルベルト空間は，サポートわれわれの問題としている測定誤差モデルにおいてはベクトルマシンをはじめとするカーネル法で幅広く使わ事後分布が閉じた形で求められないだけでなく，関数へれるようになった関数クラスであり，確率モデルであるの実質的な入力は潜在変数 z となるため，この方法は正規過程 (GP: Gaussian Process) と等価なものである

使えず、カーネル関数の基底が無限個必要となるように[6, 7]．
再生核ヒルベルト空間はガウスカーネルなど無限個の思われる．本稿で提案する手法はその点をうまく解決し∗ 産業技術総合研究所脳神経情報研究部門, 305-8568 茨城県つくば市
梅園 1-1-1 中央第 2, tel. 029-861-5549, e-mail s.akaho@aist.go.jp,
The National Institute of Advanced Industry and Scientific Technology, Central 2, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568† 統計数理研究所, 〒 190-8562 東京都立川市緑町 10-3, tel. 0505533-8500(代), e-mail iba@ism.ac.jp,
Institute of Mathematical Statistics, 10-3 Midoricho, Tachikawa,Tokyo 190-8562

張するものである．提案手法の中心となるトリックはマて、測定誤差モデルに対してカーネル回帰の考え方を拡
ルコフ連鎖モンテカルロ法を適用する際に潜在変数の次の値の候補の生成と関数のサンプリングの順序を交換することである．これによって、MCMC の各ステップでは有限次元の計算を行いつつ、全体としては無限自由度の空間をサンプルすることが可能になる．

は困難であるため，ほかの変数を固定して 1 変数ずつ更2をランダムに初期化した上で，i = 1, 2, . . . , n について新する手法を用いる．すなわちまず，z = (z1 , z2 , . . . , zn )生成モデル
入力 x, 出力 y, 潜在変数 z に関して以下の生成モデ

ルを仮定する．

p(y | z) = N [f (z), σy2 ],p(x | z) = N [0, σx2 ].

1. z を固定して f をサンプリングする
2. f と z −i を固定して zi をサンプリングする．

(2)という２ステップで交互にサンプリングを行うことを一

ただし，N [μ, σ 2 ] は平均 0, 分散 σ 2 の正規分布を表す．

つのサイクルとし，サイクルを繰り返すことによって fz の事前分布は（計算可能であれば）任意の確率分布でよいが，本稿では一様分布を仮定する．のサンプル集合を得る．

また，本稿では簡単のため上記のハイパーパラメータ

なお，以下の節で見るように，f のサンプリングはGibbs サンプラーと同様に目標分布の事後分布で生成さσx , σy は既知とする．x にのるノイズと y にのるノイせることができるが，zi の方は複雑な分布となるなのズはここでは独立であると仮定したが，相関がある場合で Metropolis-Hastings 法を用いる．

もアルゴリズムは容易に拡張できる．
一方，関数 f の事前分布はカーネル関数 k(z, z  ) をCov[f (z), f (z  )] = k(z, z  )f のサンプリングのためのトリック

z を固定して f をサンプリングする際，f は一般に共分散関数とする平均 0 の正規過程

E[f (z)] = 0,3.1

無限次元となる．単純なカーネル回帰では，入力点が固(3)定であるため Representer 定理によってその固定されたとする． カーネル関数は正定値である必要があり，例入力点だけの関数として書くことができる．ただし，測えばガウスカーネル

定誤差モデルでは入力点が潜在変数となっており，サンプリングによって時々刻々と変動するため，入力点の可k(z, z  ) = exp(−β(z − z  )2 )(4)などがある．
以上の生成モデルに基づいて n 個の入力と出力の組

D = (x1 , y1 ), (x2 , y2 ), . . . , (xn , yn ) が与えられたときにf (や (xi , yi ) に対応する潜在変数 zi ) の事後分布を求めることが問題である．なお，f, z1 , . . . , zn の事後分布は
p(f, z1 , . . . , zn | D) ∝(yi − f (zi ))2(zi − xi )212exp −(5)−−f
2σy22σx22と書ける． ただし，f  は k(z, z  ) の定める再生核ヒルベルト空間のノルムである．

3マルコフ連鎖モンテカルロ法
マルコフ連鎖モンテカルロ法 [10] は，事後分布を解析的に求めることが困難な場合にマルコフ連鎖を用いたシミュレーションによってモンテカルロサンプリングを行う手法である．マルコフ連鎖は，その定常分布が目的の事後分布となるように設計され，Metropolis-Hastings
アルゴリズムや Gibbs サンプラなどのアルゴリズムが知られている．

能性が無限に増えてしまう．例えば，f の値を，ある潜在変数 zi に対して生成した f (zi ) の値は，zi が別の値zi に変動してしまうともはや zi の関数値とみなすことはできない．
ただし，ほかの変数を固定して zi だけを変動させる
とき，それを採択するかどうかを決めるのに必要なのは後で見るように，f (zi ) と f (zi ) の二つの値だけである．

そこで，zi のほかに，zi が次に変動する候補値 zi をあらかじめある確率分布によって生成しておき，f (zi ) とf (zi ) をランダムに生成し，その後で zi を採択するかどうかを決めるというトリックを使う．これは一見 ziと f の生成する順番が逆になってように見えるが，f のサンプリングは入力点とは独立して行うことができるた
め，このようなテクニックが使えるのである．
これにより，（一般に無限自由度の）操作:「関数 f を事後分布から生成し，f (z1 ), . . .,f (zn ), f (zi ) という値を評価すること」
が，有限次元の操作:「zi の生成を f よりも先に行うことによって，f (z1 ), . . ., f (zn ), f (zi ) の事後分布から我々の問題では，D が与えられたときに潜在変数 z1 ,サンプルを生成すること」

z2 , . . ., zn と f に関してシミュレーションを行う必要があに還元できることを意味している．

るが，同時に複数の確率変数に関する乱数を生成するのzi の提案分布はエルゴード的であれば任意でよいが，本稿では xi の近くに zi があることを想定し，
1 

2p(zi ) ∝ exp − (zi − xi )2(6)に従って zi を生成する．
事後分布に従って f を生成する具体的な方法は 3.3
節で述べる．

3.2OO


(12)

成分のブロックに分割すると⎛K1 k 2
⎜ TK = ⎝ k2 k3k4⎞
K3⎟
kT4 ⎠(13)K5となる．このうち zi を zi に更新する場合に新たに計算

zi の採否プリングを行う前に決めていたが，f (zi ), f (zi ) を生成した後これを採用するかどうか決める手続きが必要となる．本稿では (6) 式のように zi を zi とは独立に生成したので Metropolis-Hastings 法による採択確率は
(yi − f (zi ))2 − (yi − f (zi ))2α = exp −2σy2が必要なのは k2 , k3 , k4 の n+1+m 個の成分だけである．

K5 は最初に一度だけ計算しておけばよく，K1 , K3 も最初に一度計算して，zi が採択された場合には，k2 , k3 , k4を zi に相当する成分 (i 行と i 列) に置き換えるだけでよい．
さて，これらの値をもとに f を実際にサンプリング(7)するのだが，そのために平均 µ, 分散共分散行列 V を求めよう．まず平均は，上記の事後分布からも導出できとなり，min(1, α) の確率で zi を採択すればよい．

るが，より計算効率がよい方法は以下の事実を利用したものである．

f のサンプリングの詳細

ここでは f のサンプリングを行う数値計算の詳細について述べる．
関数 f の統計量を計算するためには，サンプリングした f 自身を保存する必要があるが，先にも書いたよ
うに f は有限個の点でのみ評価可能な関数である．そこで本稿では，実軸上にグリッド点 u1 , u2 , . . . , um という有限個の点を追加し，これらの点で f の値をサンプ

命題 1 z を固定したときの f の事後分布の平均 (=正規分布の場合 MAP 解) は，(z1 , y1 ), . . . , (zn , yn )に対する正則化付き最小二乗カーネル回帰の解
nf (z) =i=1 αi k(zi , z) に z の各要素を代入してできるベクトルに等しく，
⎞⎛
K1
⎟
⎜
µ = ⎝ kT2 ⎠ (K1 + σy2 In )−1 y(14)K3Tリングして統計量を計算する． このように有限個グリッド点での値の評価を行っても，サンプリングしている関数全体は無限次元のままであることに注意する．
最終的に，与えられたデータ点と新規に生成した候補点，それから関数値の評価のためのグリッド点を合わせた n + 1 + m 個の点z = (z1 , . . . , zn , zi , u1 , . . . , um )T(8)における f の値が出てきた． それらをまとめてとおくと，z を固定した元での f の事後分布は

p(f | z, D) ∝11exp − 2 (y + − f )T J(y + − f ) − f T K −1 f2σy2(10)

これは Representer 定理からも理解できる．
次に分散共分散行列は (10) 式の 2 次形式を取りだして −2 倍して逆行列を取ればよく，

−11V = K −1 + 2 Jσy(15)

算結果が不安定になりやすい．そのため K の Cholesky
分解を K = LT L としたときに，

−11T
T
In+1+m + 2 LJLV =LLσy(16)と変形したものを使った方が数値的に安定であり，さらに L を K の不完全 Cholesky 分解で近似すれば計算量も減らすことができる．
結果として，平均 µ と分散共分散行列 V を計算し，N [µ, V ] に従う多次元正規乱数を発生して f のサンプという多変量正規分布となる． ただし，y = (y1 , . . . , yn )T ,で与えられる．となる． ただし通常 K は 0 に近い固有値が多いため計f = (f (z1 ), . . . , f (zn ), f (zi ), f (u1 ), . . . , f (um ))T (9)y + = (y T , 0, . . . , 0)T ,Oであり，K は z に対するグラム行列で，n, 1, m 個のK3T

zi が次にどこに動くかという候補点 zi は f のサン

3.3

J=In(11)ルとすればよい．

3.4

アルゴリズムのまとめ

以上説明してきたサンプリングアルゴリズムをまとめ2ると以下のようになる．

TargetRegression
PostMean
PostStd

1. 初期化：z1 , . . . , zn を生成する ((6) 式)．また，(13)1式のグラム行列 K のうち，K1 , K3 , K5 を計算

y2. i = 1, 2, . . . , n について以下を実行する

0する．

(a) zi を生成する ((6) 式)−1(b) グラム行列 K の k2 , k3 , k4 を計算する．

−2(c) µ, V を計算する ((14),(16) 式)．
(d) N [µ, V ] に従って f を生成する．
(e) (7) 式の確率で zi を採択し，採択されればグラム行列を更新し zi を zi で置き換える．
3. ステップ 2 のサイクルを反復し，グリッド点における f を生成する．

4−6−4−20246x

図 1: あてはめの例．σx2 = σy2 = 0.12 ．提案手法はモン
テカルロ法によるグリッド点 ({ui }mi=1 ) における f (ui )数値実験の期待値と標準偏差． 点線 (Target): データを生成

人工的に生成したデータを用いて数値実験を行った．

した f , 破線 (Regression): 通常のカーネル回帰, 実線

データはまず (−6, 6) 区間に一様分布で z を n = 20 点(PostMean): 提案手法（期待値），細い実線 (PostStd):生成し，各 zi から分散 σx2 のノイズを加えて xi を生成

エラーバー（期待値 ± 標準偏差）した．一方 yi は関数通常のカーネル回帰との比較を行った．結果は図 3 に示

f (zi ) = exp(−zi2 /5) sin(zi )(17)に分散 σy2 のノイズを加えて生成した．カーネル関数にはいずれも β = 1 のガウスカーネルを用いた．関数を評価するグリッド点は (−6, 6) の間に 50 個等間隔に置いた．

σx2 = σy2 = 0.12 のときの実行例 (図 1) を示す．提案

したとおり，ほとんどの場合でカーネル回帰の結果を上回っている．
．
次に，ノイズレベルを σx2 = σy2 = 0.32 とかなり大きくした場合にどうなるかについて同様の実験をやった結果が図 4,5,6 である．図 3 と図 6 を比べてみるとノイズレベルが高い方が提案手法による改善が大きい様子がわかる．

手法と比較する手法として，入力へのノイズを仮定しない通常のカーネル最小二乗回帰によるあてはめを行った．
提案手法の結果は繰り返しのサイクルを 100 サイクル行った単純期待値であり，burn-in 等の処理は行っていない．ノイズレベルが比較的小さいのでカーネル回帰との違いがわかりにくい (提案手法の実線とカーネル回帰の破線が重なっている部分が多い) が，この図では左端のデータがない部分でカーネル回帰のあてはまりが悪い．そこで，フィッティングを定量的に評価するために，データを生成した曲線との区間 (−6, 6) での二乗距離 (グリッド点における関数値の二乗平均) の変化をプロットしたのが図 2 である．この例では二乗距離で 25% 程度の改善が見られる．
これらのデータはあくまでワンショットのあてはめなので，乱数の種を変えて 50 回行ったときの提案手法と5おわりに測定誤差モデルにおいて，再生核ヒルベルト空間上の関数の事後分布を生成するようなマルコフ連鎖モンテカルロ法を提案した．サンプリングの順序を入れ替えるトリックによって無限次元空間のサンプリングを有限次元のサンプリングに変換することができた．
測定誤差モデルは非常に推定が困難なモデルとされている [4]．今後は理論的あるいは実験的に収束性や近似精度の分析が必要である．また，位相応答曲線の推定など実データを用いた有効性の検証も今後の課題である．
謝辞ある．

本研究は科研費 18500184 の助成を受けたもので2−20.00−10y0.02
0.01

MSE

10.03

TargetRegression
PostMean
PostStd

020406080

100−6−4−20

iteration

MSE

0.00 0.02 0.04 0.06 0.08 0.10 0.12

0.08
0.06
0.04

proposed

0.02

0.06

0.08
020

4060

80

regression
iteration

図 3: カーネル回帰と提案手法との比較 (50 回, σx2 =σy2 = 0.12 )6

図 4: あてはめの例 (σx2 = σy2 = 0.32 )カーネル回帰による結果

0.04

4x

図 2: 二乗誤差の推移 (σx2 = σy2 = 0.32 )．破線は通常の0.022

図 5: 二乗誤差の推移 (σx2 = σy2 = 0.32 )100

[8] A. J. Smola, S. V. N. Vishwanathan, and T. Hofmann. Kernel methods for missing variables. In
Proc. of 10th Int. Workshop on AI& Statistics
0.15

2005, pp. 325–334, 2005.

分析. 信学技報 PRMU2008-121, Vol. 108, No. 327,
pp. 69–74, 2008.

0.10

proposed[9] 藤木, 赤穂. 入力空間での計量に基づいた核主成分[10] 伊庭, 種村, 大森, 和合, 佐藤, 高橋. 計算統計 II. 統
0.05

計科学のフロンティア 6. 岩波書店, 2005.

0.05

0.10

0.15

regression

図 6: カーネル回帰と提案手法との比較 (50 回, σx2 =σy2 = 0.32 )参考文献[1] K. Nakae, Y. Iba, and T. Aoyagi. Statistical estimation of phase response curves and applicationin neural science. In Int. Soc. for Bayesian Analysis 9th World Meeting, 2008.
[2] 中江, 伊庭, 青柳. レプリカ交換モンテカルロ法による位相応答曲線の統計的推定. 第 11 回情報論的学習理論ワークショップ (IBIS), 2008.

[3] P. J. Bickel, C. A. J. Klaassen, Y. Ritov, and J. A.
Wellner. Eﬃcient and Adaptive Estimation forSemiparametric Models. Springer-Verlag, 1993.
[4] 甘利, 川鍋. 線形関係の推定−最小 2 乗法は最良であるのか? 応用数理, Vol. 6, No. 2, pp. 96–109,1996.
[5] S. M. Berry, R. J. Carroll, and D. Ruppert.
Bayesian Smoothing and Regression Splines for
Measurement Error Problems. J. of the American Statistical Association, Vol. 97, No. 457, pp.
160–169, 2002.
[6] 赤穂. カーネル多変量解析 — 非線形データ解析の新しい展開. 確率の情報と科学. 岩波書店, 2008.
[7] C. E. Rasmussen and C. K. I. Williams. GaussianProcesses for Machine Learning. MIT Press, 2006.