計測自動制御学会論文集
Vol.39, No.5, 503/512（2003）

変分法的ベイズ推定法に基づく正規化ガウス関数ネットワークと階層的モデル選択法
吉 本

潤一郎∗ ・石 井信∗∗ ・佐 藤 雅 昭∗∗∗Normalized Gaussian Network Based on Variational Bayes Inference and Hierarchical Model Selection
Junichiro Yoshimoto∗ , Shin Ishii∗∗ and Masa-aki Sato∗∗∗This paper presents a model selection method for normalized Gaussian network (NGnet). We introduce ahierarchical prior distribution of the model parameters and the NGnet is trained based on the variational Bayes(VB) inference. The free energy calculated in the VB inference is used as a criterion for the model selection.
In order to efficiently search for the optimal model structure, we develop a hierarchical model selection method.
The performance of our method is evaluated by using function approximation and nonlinear dynamical system
identification problems. Our method achieved better performance than existing methods.
Key Words: normalized Gaussian network, variational Bayes inference, hierarchical prior distribution, model
selection, function approximation1.

リズム 5) を用いて，高速な学習が可能であることを意味しはじめにている．しかしながら，ML 法はモデルの複雑さや推定の確正規化ガウス関数ネットワーク (normalized Gaussian

からしさを考慮しないために，モデル構造（注 1）の決定が困

network, NGnet) 1), 2) は，正規化されたガウス関数を用い難であり，また，しばしば過学習の問題を引き起こす．我々て入力空間を滑らかに分割し，分割された部分空間ごとに線は以前の研究 6) において確率的解釈に基づくユニットの生形近似を行うモデルである．空間全体をすべてのユニットを

成・削除機構を導入することによりこの問題の解決を図った

用いて近似する多層パーセプトロンのような大域モデルとが，ヒューリスティクスであるために，課題に応じてあらか異なり，NGnet は各ユニットに割り当てられた部分空間内

じめ適切なメタパラメータを設定する必要があった．

で近似を行う局所モデルであるため，効率良い学習が可能で本論文では，この問題点を解決する手法としてベイズ推定ある．局所モデルを高次元データに適用する場合には『次元法に基づく NGnet の学習法とモデル選択法を提案する．ベの呪い』が問題になるが，現実のデータは高次元空間のあるイズ推定法では，モデルパラメータに対する事前知識を事前一部分に局在しやすいため，多くの場合は局所モデルで十分布の形で導入することにより，学習データに対する各モデ分に対応可能である．実際に，NGnet は非線形力学システ

ルパラメータの確信度が事後分布として獲得される．推定後ム同定問題や強化学習における関数近似器として用いられ，には，事後分布に関するアンサンブル平均を用いて予測を行

良い結果が得られている

3), 4)

うことができる．これはアンサンブル学習の一種であるため．

NGnet は入出力変数の同時分布を近似する確率モデルとに，学習データが少ない場合でも過学習を避けることがで

1), 2)．

き，近似精度が高くなる．さらに，推定の際に計算されるモ

これは，最尤推定 (Maximum Likelihood estimation, ML)

デル周辺化尤度はモデル構造に対する尤度に相当するため法の実現法である Expectation-Maximization (EM) アルゴに，これを用いてモデル構造の良さを定量的に評価することして定式化すると，混合正規分布の特殊な場合となる

ができる．
∗

科学技術振興事業団 CREST 銅谷プロジェクト∗∗ 奈良先端科学技術大学院大学情報科学研究科∗∗∗ ATR 人間情報科学研究所∗ CREST Doya Project, Japan Science and Technology
Corporation∗∗ Graduate School of Information Science, Nara Instituteof Science and Technology∗∗∗ ATR Human Information Science Laboratories（Received July 22, 2002）（Revised February 7, 2003）

一方で，NGnet のような非線形混合モデルに対するベイズ推定は，困難な積分計算を含んでいるために近似計算が必要となる．近年，この近似手法として変分法的ベイズ推
定 (Variational Bayes inference, VB) 法が提案された 7)〜9) ．

VB 法では，ある試験事後分布を用いて真の事後分布の近似を行う．この近似は，両分布間の Kullback-Leibler (KL)
（注 1） NGnet ではユニット数がモデル構造に対応する．

c 2002 SICETR 0005/103/3905–0503 °504

T. SICE Vol.39 No.5 May 2003

divergence によって定義される自由エネルギーの最大化問

Pは混合比であり，gi ≥ 0 (i = 1, · · · , M )，かつ， M
i=1 gi = 1

題を解くことによって実現される．推定後の自由エネルギーを満たす．NN (x |µi , Si ) は N 次元中心ベクトル µi と N ×N

は (対数) モデル周辺化尤度の下界になっているので，最適次逆共分散行列 Si をパラメータとして持つガウス関数をモデル構造を決定するための評価基準を与える．また，試験表している (付録 A. 1 を参照)．Wi ≡ (wi1 , · · · , wiD )0 は事後分布における未知変量の部分独立性を仮定すると，VB

D × (N + 1) 次線形回帰行列を表すパラメータであり，wih

法は ML 法における EM アルゴリズムと類似した効率良いは出力 yh に対応する N + 1 次元線形係数ベクトルである．

NGnet は入力と出力の組 (x, y) を確率事象とする確率モ

反復アルゴリズムとして実装することができる．

VB 法における反復アルゴリズムは自然勾配法 10) の一種

デルとして定式化できる 1), 2) ．各事象は 1 つのユニット i かである 9) ために高速な収束が期待される一方で，アルゴリ

ら確率的に生成されると仮定し，生成したユニットインデク

ズムの初期条件に依存して局所最適解に収束する．この局所

ス i を隠れ変数とみなす．このとき，(x, y, i) の組に対する最適解の良さは，さらにモデル構造にも依存する．したがっ確率モデルを以下で与える．

て，悪い局所解を避けながら良い結果を得るためのモデル構
造探索手法は実用上重要な問題であり，本研究では階層モデ
ル選択法を導入することによりこの問題の解決を試みる．
本研究と独立・並行して，Ueda らは一般の混合モデルに対して，VB 法に基づくモデル選択法を提案している 11) ．
ここでも，その応用例として NGnet が取り上げられているが，本研究と以下の点で異なる．第一にモデルパラメータに関する事前分布が異なる点である．Ueda らの研究では，

Automatic Relevance Determination (ARD) 12) と呼ばれる事前分布のパイパーパラメータのみに階層的な事前分布を与えていたが，本研究では入出力分散のハイパーパラメー
タにも階層的な事前分布を与える．これにより，事前分布に

P (x, y, i|θ) = gi NN (x |µi , Si ) ND (y |Wi x̃, Bi )ここで，Bi ≡ diag (βi1 · · · , βiD ) は出力逆分散を表すパラメータである．θ ≡ {(gi , θi )|i = 1, · · · , M } は NGnet のモ
デルパラメータの集合であり，θi ≡ {µi , Si , Wi , Bi } である．
この確率分布から，入力 x が与えられた時の出力 y の期待
R
dyP (y|x, θ)y を求めると，(1) 式で定義され

値 Eθ [y|x] ≡た NGnet の出力と一致する．すなわち，確率分布 (2) 式はNGnet の確率モデルを定義している．
(2) 式より，観測変数の組 (x, y) に対する同時分布は以下の混合モデルとして定式化できる．
X
P (x, y|θ) =P (x, y, z|θ){z}

よる推定のバイアスを極力抑えることができる．第二にモデル構造探索アルゴリズムが異なる点である．Ueda らの手(2)

P (x, y, z|θ) = exp

法は，ユニットを局所的に分割・併合することによりユニッ

"M
X

(3a)
#

zi log P (x, y, i|θ)

(3b)

i=1

ト配置の不均衡性を解消する局所的なアプローチである．一

ここで，z ≡ (z1 , · · · , zM )0 は M 項変数であり，観測変数方で，本研究で提案する手法は，データ空間全体を部分空間に階層的に分割し，階層構造を訪問しながら各部分空間で最(x, y) がユニット i から生成されたことを zi = 1，zj = 0
P(j 6= i) によって表す．定義より M
i=1 zi = 1 である．(3) 式

適化を行う手続きであるので，大域的なデータ構造も考慮しによる定式化では，隠れ変数は i から z に変換されている．

たものになっている．

3.

本論文は以下のように構成される．第 2 節で NGnet とそ学習アルゴリズムの確率モデルについて概説する．第 3 節では，階層的事前

3. 1

分布を導入し，VB 法に基づく NGnet の学習法を導出する．

T 個のデータセット (X, Y ) ≡ {(x(t), y(t))|t = 1, · · · , T }第 4 節では，VB 法に基づくモデル構造選択法を提案し，そが観測された時，ベイズ推定の目的は未知変量に関する

ベイズ推定法とモデル構造評価の実装法について議論する．第 5 節では，提案手法を関数近事後分布 Ppost (Z, θ|X, Y ) を求めることである．ここで，似問題と未知の非線形力学システムの同定問題に応用する

Z ≡ {z(t)|t = 1, · · · , T } は学習データの系列に対応する隠ことにより，その性能を示し，第 6 節で本論文をまとめる．

れ変数の系列である．ベイズの定理から，未知変量に関する

2.

事後分布は以下で与えられる．

NGnet と確率モデル

Ppost (Z, θ|X, Y ) =N 次元入力ベクトル x ≡ (x1 , · · · , xN )0 を D 次元出力ベ
クトル y ≡ (y1 , · · · , yD )0 に変換する NGnet は以下で定義される

1), 2)QT

t=1 P (x(t), y(t), z(t)|θ) は完全

データセット (X, Y, Z) に対するモデルパラメータ θ の尤度．

y=

ここで，P (X, Y, Z|θ) ≡P (X, Y, Z|θ)P0 (θ)
P (X, Y )

M
X
i=1

!
gi NN (x |µi , Si )PM
j=1 gj NN (x |µj , Sj )

Wi x̃

(1)

である．P0 (θ) はモデルパラメータ θ の事前分布である．正RP{Z} dθP (X, Y, Z|θ)P0 (θ) はモデル

規化項 P (X, Y ) ≡周辺化尤度と呼ばれている．

0

ここで，x̃ ≡ (x0 , 1) であり，プライム記号 (0 ) は転置を表

ベイズ推定において計算されるモデル周辺化尤度は，最適している．M は NGnet を構成するユニット数であり，i ∈なモデル構造，すなわち，NGnet のユニット数 M を決定

0

するための定量的な評価基準を与える．モデル構造 M への{1, · · · , M } はユニットインデクスである．g ≡ (g1 , · · · , gM )

計測自動制御学会論文集 第 39 巻 第 5 号

2003 年 5 月

505

依存性を明確にするためにモデル周辺化尤度を P (X, Y |M )パラメータは定数である．この定式化では，ハイパーパラと記述すると，これは観測データセットに対するモデル構造

メータ ξ も未知変量であるので，ベイズ推定は事後分布

M の尤度を意味している．したがって，P (X, Y |M ) が最大

Ppost (Z, θ, ξ|X, Y ) を求めるものとして拡張される．このよとなるモデル構造 M を選択することはモデル構造に関する

うな階層的事前分布を導入する利点については，3. 4 節で述最尤推定に対応する．

べる．

しかしながら，NGnet のような非線形混合モデルに対し

3. 3

変分法的ベイズ推定法て，積分計算を含むモデル周辺化尤度を解析的に求めるこVB 法 7)〜9) では，ある試験事後分布 Q(Z, θ, ξ) を用いてとは困難であり，近似手法が必要である．以下では，自然共真の事後分布 Ppost (Z, θ, ξ|X, Y ) の近似を行う．この近似

役分布を用いた階層的な事前分布を導入し，VB 法に基づいは，以下で定義される自由エネルギーの最大化問題を解くこて，パラメータ事後分布とモデル周辺化尤度を近似するためとによって実現される．
flfi
P (X, Y, Z|θ)P0 (θ|ξ)P0 (ξ)
F [Q] ≡ log
Q(Z, θ, ξ)Qのアルゴリズムを導出する．

3. 2

階層事前分布

本研究では，モデルパラメータの事前分布として以下で定義される自然共役分布を与える．

P0 (θ|ξ) = P0 (g)MY

下で定義される．

P0 (µi |Si )P0 (Si |σi )P0 (g) = DM (g|γ0 )

(4a)
(4b)

P0 (µi |Si ) = NN (µi |m0i , γ0i Si )(4c)

P0 (Si |σi ) = WN (Si |γs0 , γs0 σi IN )(4d)

P0 (Wi |Bi , Υi ) =DY

NN +1 (wij |0, βij Υi )

(4e)

D
Y

G (βij |γβ0 /2, γβ0 ρij /2)(4f)

ここで，DM (·|·)，WN (·|·, ·)，および，G (·|·, ·) は，それぞれ，M 次元 Dirichlet 分布，N 次元 Wishart 分布，および，ガンマ分布を表記している (付録 A を参照)．IN は

N × N 次単位行列を表記している．γ0 ≡ (γ01 , · · · , γ0M )0であり，γs0 と γβ0 はスカラーである．m0i は N 次元ベ
クトルであり，事前分布における µi の期待値に対応す`´
る．Υi ≡ diag υi1 , · · · , υi(N +1) は，事前分布におけるWi の逆分散を制御するパイパーパラメータである．σi とRi ≡ diag (ρi1 , · · · , ρiD ) は，それぞれ，事前分布におけるSi と Bi の期待値の逆数を表現するハイパーパラメータである．Ueda らの研究 11) では，以上のハイパーパラメータのうち Υi のみを未知変量として扱い，階層的事前分布を与えているが，本研究では，ξ ≡ {σi , Υi , Ri |i = 1, · · · , M } を未知変量として扱い，以下の階層事前分布を与える．

P0 (ξ) =P0 (σi )P0 (Υi )P0 (Ri )

(5a)

i=1
−1
P0 (σi ) = G(σi |γσ0 /2, γσ0 τσ0/2)

P0 (Υi ) =

N
+1
Y

(5b)

D
Y

dθdξQ(Z, θ, ξ)f (Z, θ, ξ){Z}

KL (Q kPpost ) は，試験事後分布 Q と真の事後分布 Ppost のKL divergence であり，以下で定義される．
flfi
Q(Z, θ)KL (Q kPpost ) ≡ log
Ppost (Z, θ|X, Y ) Q

第 1 項の log P (X, Y ) が Q に依存しないことに注意すれば，
自由エネルギー F [Q] を Q に関して最大化すると，試験事
後分布 Q は真の事後分布 Ppost に等しくなり，最大化された後の自由エネルギーは最適モデル構造選択基準である（対数）モデル周辺尤度と等しくなることが分かる．
本研究では，隠れ変数 Z ，モデルパラメータ θ，およびハ
イパーパラメータ ξ が互いに独立となるような試験事後分
布 Q を用意し，それによって事後分布 Ppost (Z, θ, ξ|X, Y )の近似を行う．すなわち，Q(Z, θ, ξ) = QZ (Z)Qθ (θ)Qξ (ξ)と因子分解できるものとする．この時，(6) 式の自由エネル
ギーは以下のように変形できる．

F [QZ , Qθ , Qξ ] = L − (H θ + H ξ )
fl fl
fifi
P (X, Y, Z|θ)L≡
log
QZ (Z)
θ Z*fifl +Qθ (θ)Hθ ≡log
P0 (θ|ξ) ξ
θ
fi
fl
Q
(ξ)
ξ
H ξ ≡ log
P0 (ξ) ξ

(7a)
(7b)
(7c)
(7d)

ここで，
−1
G(υij |γυ0 /2, γυ0 τυ0
/2)

(5c)

j=1

P0 (Ri ) =XZ

布として一致する時に，最小値 0 になる．また，(6) 式右辺

j=1

M
Y

hf (Z, θ, ξ)iQ ≡KL divergence は，Q(Z, θ) と Ppost (Z, θ|X, Y ) とが確率分

j=1

P0 (Bi |Ri ) =(6)

ここで，h·iQ は試験事後分布 Q に関する期待値であり，以i=1×P0 (Wi |Bi , Υi )P0 (Bi |Ri )= log P (X, Y ) − KL (Q kPpost )hf (Z)iZ ≡

X
{Z}

−1
G(ρij |γρ0 /2, γρ0 τρ0
/2)

(5d)

j=1

以上の事前分布において，添字 “0 ” を持つ全てのハイパー

QZ (Z)f (Z),Z
hf (θ)iθ ≡ dθQθ (θ)f (θ),
の表記法を用いた．

Z
hf (ξ)iξ ≡ dξQξ (ξ)f (ξ)

506

T. SICE Vol.39 No.5 May 2003

(7) 式において，L は期待対数尤度に対応し，観測データ
θに対するモデルの適合度を表す項である．一方，H と H

ξして扱い，階層的な事前分布を与えたが，これには VB-EMH
アルゴリズムの実装の観点から以下のような利点がある．

は，それぞれ，モデルパラメータとハイパーパラメータに関（ 1 ） 各ユニットの線形回帰行列 Wi は D × (N + 1) 個のするエントロピーに対応し，これらはモデルの複雑さに対す

自由パラメータから構成される．一方で，ユニットによっ8), 9)．したがって，VB 法はモデルのては出力回帰に無関係な入力次元が存在する場合がある．

複雑さに対する正則化を導入した期待対数尤度の最大化手

この場合に少ない観測データから Wi を推定すると過剰な法とみなすことができる．モデルの複雑さに対する正則化項

自由度のために過学習を招きやすい．これを防ぐために，の効果は，事前分布における定数ハイパーパラメータの影響

Mackey らによって提案された ARD 法 12) に基づく階層を大きく受ける．しかし，データの分布に関する事前知識が事前分布を (4e) 式と (5c) 式によって導入している．このない場合，これらのハイパーパラメータを直接設定し，正則

階層事前分布は，事前分布に関する Wi の期待値を 0 と化項の効果を制御することは困難な場合がある．そこで本研し，その逆分散を Υi に比例させて与えるものである．ベ

究では，より積極的にモデルの複雑さに対する正則化を制

イズ推定を用いて Υi を推定すると，出力回帰に無関係な御するために，確信度 κ (κ > 0) を導入する．これにより，

入力 xj に対応する逆分散 υij は非常に大きくなる．これ(7) 式の自由エネルギーは以下で修正される．

により，υij に対応する回帰パラメータは実質的に定数とるペナルティ項となるF κ [Qz , Qθ , Qξ ] ≡ κL − (H θ + H ξ )(8)

みなせるので，Wi の自由パラメータ数が自動的に決定できることになる．

(8) 式はデータを κ 回重複して観測したことに対応し，VB（ 2 ） (4d) 式における γs0 は，入力逆共分散行列 Si に関法における観測データに対する信頼度が通常に比べて κ 倍する事前知識がデータ何個分に相当するかを決める定数になっていることを意味している．

ハイパーパラメータである．一方，入力逆共分散行列を推自由エネルギー (8) 式の最大化は，以下の 3 ステップを反復手続きによって行う（実装法は付録 B を参照）．

1. VB-E ステップ:とに対応して，Wishart 分布は γs0 ≤ N − 1 の時に発散

Qz ← argmax F κ [Qz , Qθ , Qξ ]Qz

定するためには少なくとも N 個のデータが必要になるこする．そのために，γs0 は (N − 1) より大きくなければな

2. VB-M ステップ:Qθ ← argmax F [Qz , Qθ , Qξ ]らない．特に入力次元 N が大きい場合，この制限はパラ

3. VB-H ステップ:Qξ ← argmax F κ [Qz , Qθ , Qξ ]メータ事後分布の推定に非常に強いバイアスを加えること

κQθを意味する．この困難を避けるために，我々は Si の事前

Qξ

上記の反復手続きを ML 法における EM アルゴリズム 5), 13)（注 2）との類似性から VB-EMH アルゴリズム

分布としてその期待値が σi−1 IN となる Wishart 分布を与と呼ぶ．VB-

えている．この事前分布は Si に関する事後分布を推定す

EMH アルゴリズムの各ステップで自由エネルギーは単調増る際に hSi iθ が縮退するのを防ぐ正則化の役割を果たし，加するので，自由エネルギーは局所最大値に大域的に収束す

数値計算上安定した結果を得ることができる（付録 B にる 9) ．

おける (B. 1) 式と (B. 2) 式を参照）．さらに，σi を未知学習後，NGnet の出力 ŷ は予測分布に基づいて以下で与えることができる．
Zŷ =dydθPpost (θ|X, Y )P (y|x, θ)yZ≈
dydθQθ (θ)P (y|x, θ)y = hEθ [y|x]iθよるバイアスが観測データに適応して自動的に決定される．これによって推定における事前分布の依存性を極力抑えることが期待できる．

(9)

ML 法や最大事後確率推定（Maximum A Posteriori estimation, MAP）法では，最適なモデルパラメータを θ∗ として点推定し，NGnet の出力は Eθ∗ [y|x] となる．一方で，VB 法では (9) 式のように各モデルパラメータ θ に対するEθ [y|x] が事後分布で重み付き平均された形が NGnet の出力となる．これはアンサンブル学習の一種であり，ML 法や

MAP 法に比べて過学習を防ぐ効果が期待でき，かつ，近似精度が高くなる．

3. 4

変量として扱いベイズ推定を行うことによって，正則化に

階層事前分布の利点

本研究では，事前分布ハイパーパラメータ ξ を未知変量と（注 2） 通常，VB-EM アルゴリズムと呼ばれるが，ハイパーパ
ラメータに関する推定ステップを含んでいるために，ここではVB-EMH アルゴリズムと呼んでいる．

（ 3 ） 上記と同様の理由により，出力誤差逆分散 Bi の推定に関わるハイパーパラメータ ρij に対しても階層事前分布を導入し，ベイズ推定によって事前分布によるバイアス
を観測データに基づき自動的に決定する．
上記の第 (1) 項については Ueda らの研究 11) でも導入されているが，第 (2) 項と第 (3) 項については本研究で初めて導入されたものである．

4.

階層的モデル選択法

VB-EMH アルゴリズムは局所的な最適解に収束するので，NGnet のような非線形モデルでは試験事後分布の初期値によって収束先が異なる．この局所最適解の良さは，さらにモデル構造にも依存する．そのため，悪い局所最適解を避けながら良いモデル構造を選択するための手法は実用上重要である．本節では，データセットの大域的な分布を考慮に入れ

計測自動制御学会論文集 第 39 巻 第 5 号

2003 年 5 月

507

た階層的モデル選択法を提案する．この手法は，NGnet のみならず一般の混合モデルに対しても適用可能である．

4. 1

実装法

本手法は 2 つのフェーズから構成される．第 1 フェーズ
（トップダウンフェーズ）では，学習データセットを階層的に分割し，2 分木を構成する．根ノードでは M = 2 個のユニットを持つモデルを用意し，VB-EMH アルゴリズムを用いて学習データセット S ≡ (X, Y ) に対する学習を行う．

VB-E ステップにおいてユニットインデクス i に関する事後分布を得ることができるので，その MAP 基準に従って Sを 2 つのサブセット S1 と S2 に分割する．すなわち，Si は，

i = argmaxj Qz (zj (t) = 1) を満たすデータ (x(t), y(t)) ∈ Sの集合である．分割されたデータサブセット S1 ，S2 のぞれ
ぞれに対して上記の手続きを再帰的に行うことにより，データについての 2 分木が構成できる．なお，この再帰手続き
は，与えられたデータサブセットについて M = 1 個のユ

Fig. 1 Binary data tree produced in top-down phase

ニットを持つモデルが M = 2 個のユニットを持つモデルよりも局所自由エネルギー（そのデータサブセットに対してのみの自由エネルギー）が大きくなったところで停止する．図

1 は 2 次元混合正規分布の推定に階層的モデル選択法を適用した時のトップダウンフェーズの様子を示したものである．
ここで，“4” 印と “x” 印は各階層で分割されたデータサブセットを表しており，楕円は推定に用いられた各ユニットの正規分布を表している．
第 2 フェーズ（ボトムアップフェーズ）では，トップダウンフェーズで構成された 2 分木の各ノードを下から上へと深さ優先で訪問しながら，ユニットの統合を試みる．葉ノー
ドでは，トップダウンフェーズで用いられた M = 1 個のユニットを持つモデルを保持しておく．上位ノードでは，子ノードが持つ全てのユニットを統合したモデルを用意する．
ユニット統合によって生成されたモデルは，このノードに割り当てられたデータサブセットを使って VB-EMH アルゴリ
ズムによる再学習を行う．次にこのデータサブセットに対する最適モデル構造を選択するために混合比 hgi iθ が小さいユ
ニットについて順次削除操作を行いながら，局所自由エネル
ギー最大化基準により局所最適モデル構造を決定する．この手続きは下位ノードから上位ノードに向かって順番に行われる．以上の手続きを下位ノードから根ノードに到達するまで深さ優先で行うことによって，全データセットについての最適なモデル構造を持つモデルが構成される．
以上で説明した階層的モデル選択法は図 2 で示される再帰的手続きを用いて実装することができる．ここで，ステップ(1)-(7) がトップダウンフェーズに対応し，ステップ (8)-(10)がボトムアップフェーズに対応する．

4. 2

階層的モデル選択法の特徴

混合モデルでは，一般にモデルが複雑になる（混合数 M¶³

Function HMS(Dataset S)
Output: Posterior Q∗θ（ 1 ） M = 1 のモデルを用意し，VB-EMH アルゴリズムによって S に対する学習を行う．これによって得られる自由エネルギーを F1λ とし，θ に関する試験事後分布を Q̄θとする．
（ 2 ） M = 2 のモデルを用意し，VB-EMH アルゴリズムによって S に対する学習を行う．これによって得られる自由エネルギーを F2λ とし，隠れ変数 z(t) に関する試験事
後分布を Qz とする．
（ 3 ） も し ，F1λ > F2λ な ら ば ，こ の 関 数 を 終 了 し ，Q∗θ = Q̄θ を出力として返す．そうでなければ以下の処理を行う．
（ 4 ） S1 , S2 を空集合として用意する．
（ 5 ） t = 1, · · · , T に対して以下の処理を行う．
（ a ） もし Qz (z1 (t) = 1) ≥ 0.5 ならば，S1 にデー
タ (x(t), y(t)) を加える．そうでなければ，S2 にデータ
(x(t), y(t)) を加える．
（ 6 ） 関数 HMS(S1 ) を適用し，それによって得られる試
験事後分布を Qθ,1 とする．
（ 7 ） 関数 HMS(S2 ) を適用し，それによって得られる試
験事後分布を Qθ,2 とする．
（ 8 ） θi に関しては Qθ,1 と Qθ,2 の直積として与えられ，g に関しては Qθ (g) = P0 (g) で与えられる試験事後分布
Qθ を用意し，データセット S に対する学習を行う．これによって得られる自由エネルギーを F λ とし，削除候補となるユニットインデクスを i− = argmini hgi iθ とする．
（ 9 ） Qθ からユニット i− を削除した試験事後分布 Q−
θ
を用意し，データセット S に対する学習を行う．これによって得られる自由エネルギーを F̂ λ とする．
（10） もし，F λ ≥ F̂ λ ならばこの関数を終了し，Q∗θ = Qθを出力として返す．そうでなければ，Qθ = Q−θ とし，削除候補となるユニットインデクスを i− = argmini hgi iθ
とする．その後，ステップ (9) を行う．

µ´

が大きくなる）にしたがって局所最適解の数が増加する．また，悪い局所最適解では，あるデータ領域を複雑過ぎるモデルを用いて説明している一方で，別のあるデータ領域では簡Fig. 2 Procedure for hierarchical model selection

T. SICE Vol.39 No.5 May 2003
Correlation between nMSE & FE

果たす．トップダウンフェーズでは，高々M = 2 のモデルを用いて推定するために，局所最適解は少ない．このため，Correlation between NU & FE

-3200

-3200

-3400

-3400

-3600
-3800
-4000

0.1

Correlation between NU & nMSE

0.09
0.08
0.07

-3600

nMSE

こる．トップダウンフェーズはこうした状況を避ける役割をFree Energy (FE)単過ぎるモデルを用いて説明しようすることがしばしば起Free Energy (FE)

508

-3800

0.06
0.05
0.04

-4000

0.03

-4200

-4200

-4400
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

-4400

0.02
5

10

15

20

25

30

35

40

45

50

Number of Units (NU)

nMSE

0.01

5

10

15

20

25

30

35

40

45

50

Number of Units (NU)

初期値に対する依存性が減少し，安定した解が得られることが期待できる．また，推定後の事後分布を用いてデータ分割を行うためにデータとモデルの双方に適合した形でデー

Fig. 3 Correlation diagrams among generalization error(nM SE), free energy (FE) and the number of units(NU).

タ領域分割を行うことができる．さらに，各データサブセッ
トに対しても自由エネルギーが計算できるために領域分割処理に必要な計算量が大きく増加する．一方，ユニット併を進めていくことが適切かどうかを定量的に評価すること

合操作に対応する本手法のモデル探索法は，ボトムアップができる．

フェーズのユニット削除操作である．削除を試みるユニッ

ボトムアップフェーズは，分割されたデータ領域について

トの選択基準は混合比 hgi iθ で与えられており，M 通りの最適化を行う役割を果たす．トップダウンフェーズにおけの順位付けのみを行えばよい．したがって，ユニット数

る前処理によって，あるノードに割り当てられたデータ領域

M が大きくなってもこの処理に必要な計算量はそれほどはその 2 つの子ノードに割り当てられたデータ領域の和集増加しない．

合として表現できる．このため，あらかじめ子ノードが担当• 最適モデル構造が M ∗ とすると，本手法のデータサブするデータ領域に対して最適な分布が分かっていると仮定す

セットは最低でも 2M ∗ − 1 通り生成され，それぞれに対ると，親ノードでは 2 つの分布の混合分布が準最適な分布

して VB-EMH アルゴリズムを適用する必要がある．これの 1 つであると考えられる．この考察に基づいて，本手法が本手法の最大の欠点である．しかしながら，トップダウ

は各ノードが担当するデータサブセットに対する最適分布

ンフェーズでは高々2 個のユニットを持つモデルの学習をを葉ノードから逐次的に推定し，親ノードの学習を行う際行うだけであり，下位ノードに進むにつれて学習データもには，2 つの子ノードの推定結果を混合させた分布を初期値少なくなるので，それほど計算量は多くならない．ボトムとして用いている．これにより，良い初期値から学習を行う

アップフェーズでは，多くのユニットを持つモデルに対すことができるので，悪い局所解に陥ることを防ぐとともに，る学習も必要となるが，各ユニットは下位層において局所

VB-EMH アルゴリズムの反復回数を低く抑えることができ的な最適化が行われているため，収束までの EHM アル

る．一方，子ノードでは親ノードの部分領域に関する最適化

ゴリズムの反復回数はそれほど多くならない．この点は，

しか行わないために，2 つの子ノードの領域境界付近の推定

5. 1 節の実験結果によって示す．

が悪くなる場合がある．この状況は領域境界近傍で本来は局

5.

所的に 1 つのクラスタを形成しているデータがトップダウンフェーズで 2 分割されてしまった時によく起こる．このよ
うに過剰に分割されたデータ領域を担当するユニットは他のユニットに比べて混合比が小さくなると考えられるので，混合比の小さなユニットの削除を試みる．これによって，局所的な領域分割による副作用が解消されやすいと考えられる．

5. 1

実験

関数近似問題

最初の実験では，以下の関数 14) を用いて本手法の関数近似能力とモデル選択能力を評価した．
n
o
2
2
2
2
y = max e−10x1 , e−50x2 , 1.25e−5(x1 +x2 )Ueda らのユニットの分割・併合 (Split and Merge, SM)

定義域 −1 ≤ x1 , x2 ≤ 1 から入力変数 x ≡ (x1 , x2 ) を一様

法 11) も悪い局所最適解を避けながらモデル構造探索を行うにサンプルして得られた 500 個の入出力変数の組 (x, y) を手法であるが，それと比較して本手法の長短は以下のように学習データとする．ただし，各出力変数には平均 0，標準偏まとめられる．

差 0.1 の大きなガウス雑音を付加する．

• SM 法におけるユニットの分割は本手法のトップダウン

まず，混合数 M が 5 から 50 まで 46 種類のモデル構造をフェーズで行われる処理に対応する．SM 法ではデータに用意し，各モデル構造に対して VB-EMH アルゴリズムの初

対する局所的な適合度にしたがって分割を試みるユニット

期条件を変えて 100 回学習を行った．ここで，確信度 κ = 7が選択される．一方で，本手法はデータ空間全体から逐を用いている．図 3 は，モデル構造と学習後の汎化誤差，お次・階層的な分割を行っているので，大域的な分布を考慮よび，自由エネルギーの相関図を示したものである．汎化誤した分割方法となっている．

差は，入力変数をその定義域から 41 × 41 個の点として格子• SM 法において併合を試みるユニットの選択基準は事上に等間隔でとり，そこでの出力変数に関する真値と近似値

後分布 Qz (zi (t)) のユニットインデスク i に関する相関にとの平均自乗誤差を真値の分散で正規化した値（正規化自乗

よって与えられている．この計算量は O(M C2 × T ) とな

誤差, nM SE ）によって評価されている．この図から，汎化

り，かつ，M C2 通りの組合せに対して順位付けが必要と誤差の小ささと自由エネルギーの大きさの間に強い相関がなる．そのため M や T が大きくなるにしたがって，このあることが分かる．

計測自動制御学会論文集 第 39 巻 第 5 号

2003 年 5 月

Table 1 nM SE of each method in “pumadyn8” tasks

1
1.5
0.5

VB-NGnet
OEM-NGnet
MEnet
HMEnet

1
0

0.5
0

-0.5

1

509

nm256
0.12560.21040.12700.1893nh256 nm1024 nh1024
0.4901 0.0517 0.3503
0.5152 0.0518 0.3896
0.5370 0.0461 0.3603
0.5375 0.0490 0.3639

1
-1
-1

-0.5

0

0.5

1

(b) Receptive Field

x3

0
-1 -1
(a) approximated function(nMSE=0.01439)

Fig. 4 (a) Function shape approximated by the NGnet; (b)receptive fields of units in the trained NGnet.

50

50

40

40

30

30

x3

0

20
10
0
50

20

0

x2

次に階層的モデル選択法を用いてモデル選択を行った．そ

20
10
0
50

-50 -20x2

x1(a) True System

20

0

0

0
-50 -20x1(b) Reconstructed Systemの結果，混合数 M = 18 のモデル構造が選択され，汎化
誤差は nM SE = 0.0144 となった．収束するまでに要する

VB-EMH アルゴリズムの反復回数は 4807 回であった．図
4 は学習後の NGnet が出力する関数形と各ユニットの受容野を示したものである．ユニットが定義域内に均等に配置され，その結果，精度の高い関数近似が実現されている様子が分かる．
同じ実験条件において，我々が以前に提案したオンライ
ン EM 法によって得られた結果は nM SE = 0.0189 であった 6) ．VB 法と階層的モデル選択法を用いることにより近似精度は約 24% 改善された．また，モデル選択能力を比較するために SM 法によるモデル選択アルゴリズム 11) を実装した．ここで，初期モデル構造は階層的モデル選択法とFig. 5 (a) True Lorenz attractor; (b) an attractor reconstructed by the trained NGnet

上の比較から，本手法は既存の手法に比べて関数近似能力に優れた効率良いモデル選択を実現しているといえる．

5. 2

非線形力学システムの同定問題

未知の非線形力学システムの同定問題に本手法を適用した．本実験では，以下の微分方程式で定義される Lorenz シ
ステムを取り上げる．
0
1 0
1
ẋ1
a(x2 − x1 )B
C BCB ẋ2 C = B −x2 + (b − x3 )x1 C
@
A @
A
ẋ3
−cx3 + x1 x2

同じ M = 2 とした．その結果，最善の場合でも汎化誤差

ここで，x ≡ (x1 , x2 , x3 ) はシステムの状態変数である．ẋ

は nM SE = 0.0168 であり，この最善の結果を得るためには状態変数 x の時間微分を表している．システムパラメー

要する VB-EMH アルゴリズムの反復回数は 4317 回であっタは a = 10.0, b = 28.0, c = 8/3 というしばしば用いられた．階層的モデル選択法は SM 法に比べて若干多くの計算る値に設定する．図 5(a) は Lorenz システムのアトラクタ

回数を必要とするが，精度の高い関数近似を実現できているを示したものである．
アトラクタ上で時間間隔 ∆τ = 0.01 でサンプリングされといえる．
次に，関数近似問題のベンチマーク集である DELVE デー

た T = 5000 個の時系列 X ≡ {(x(k∆τ )|k = 1, · · · , T )} をの中から “pumadyn8” という課題を用いて，高観測時系列とする．観測時系列 X に対する離散化ベクトル

次元データに対する性能を評価した．この課題は，ロボッ場は以下で定義される．
x((k + 1)∆t) − x(k∆t)
V (x(k∆t)) =∆t
NGnet は x から V (x) への写像を近似するように学習を行

タセット

15)

トアームの 8 次元の状態変数から残りの 1 次元の状態変数を推定するものである．データセットはデータ数とノイズの大きさに応じていくつかの条件に分かれているが，この中から “nm256”，“nh256”，“nm1024”，“nh1024” を実験に用いた．ここで，“nm” は中程度のノイズを，“nh” は大きなノイズを含んでいることを意味している．“256” とう．学習後には，任意の初期状態 x̂(0) から離散化された時系列を以下の式を用いて自動的に生成することができる．

x̂(t + ∆t) = x̂(t) + V̂ (x̂(t))∆t“1024” は学習データセットのデータ数である．各データセッ

ここで，V̂ (x) は NGnet によって近似された離散化ベクト

トに対して本手法（VB-NGnet）を適用した時の nM SE をル場である．

表 1 の上段に示す．また，オンライン EM アルゴリズム
6)

図 5(b) は NGnet によって再構成されたアトラクタを示や，NGnet と類似した混合回帰モデル

している．Lorenz システムが持つ奇妙なアトラクタ構造をである Mixture of Experts network（MEnet） 16), 17) ，お良く再現できていることが分かる．ここで，アトラクタ上に（OEM-NGnet）

17), 18)

による結果も表 1

おける離散化ベクトル場の汎化誤差は nM SE = 0.007% でに示す．nm1024 を除く全ての条件において，本手法は他のあった．図 6 は 100 個の異なる初期状態から NGnet によっ手法より優れており，特に，データ数が少なくノイズが大きて予測された時系列と真の時系列との誤差の平均を示してい場合に，本手法が良い性能を示していることが分かる．以いる．平均して約 t = 2 まで正確な予測を行っていることがよび，階層的 MEnet（HMEnet）

510

T. SICE Vol.39 No.5 May 2003
1.6

Prediction Error

1.4
1.2
1
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2

2.5 3 3.5 4 4.5 5

Time
Fig. 6 Average prediction error of the trained NGnet

分かる．我々は以前オンライン EM アルゴリズムを用いて同じ実験を行っている 3) ．この時，ほぼ同等の結果を得るた
めには T = 100, 000 個の観測系列が必要であった．20 分の

1 のデータ量から同程度の結果を得ていることから，提案手法が高い汎化性能を持っていることが分かる．

6.

おわりに

本論文では，正規化ガウス関数ネットワークに対する新しいモデル選択法を提案した．本手法では，モデルパラメータについての階層的事前分布が導入され，変分法的ベイズ推定法に基づいて学習が行われた．この時に計算される自由エネルギーは最適モデル構造を決定するための評価基準として用いられた．また，効率良いモデル構造選択を実現するための階層的モデル選択手法が導入された．本手法は関数近似問題と非線形力学システムの同定問題に適用された．
その結果として，本手法が近似精度が良く，汎化能力に優れていることが示された．本手法はバッチ学習のためのモデル
選択手法であるが，オンライン学習に適用できる手法へ拡張することが今後の課題である．

謝辞

mum likelihood from incomplete data via the EM algorithm, Journal of Royal Statistical Society B, 39-1, 1/38
(1997)
6）M. Sato and S. Ishii: On-line EM algorithm for the normalized Gaussian network, Neural Computation, 12-2,
407/432 (2000)
7）S. Waterhouse, D. Mackay and T. Robinson: Bayesian
methods for mixture of experts, Advances in Neural Information Processing Systems 8 (eds. D. S. Touretzky and
M. C. Mozer and M. E. Hasselmo), MIT Press, 351/357(1996)
8）H. Attias: A variational Bayesian framework for graphicalmodels, Advances in Neural Information Processing Systems 12 (eds. S. A. Solla, T. K. Leen and K. R. Müller),
206–212, MIT Press, 206/212 (2000)
9）M. Sato: Online model selection based on the variational
Bayes, Neural Computation, 13-7, 1649/1681 (2001)10）S. Amari: Natural gradient works efficiently in learning,
Neural Computation, 10-2, 251/276 (1998)11）N. Ueda and Z. Ghahramani: Bayesian model search for
mixture models based on optimizing variational bounds,
Neural Networks, 15, 1223/1241 (2002)12）D. J. C. Mackay: Bayesian non-linear modelling forthe prediction competition, ASHRAE Transactions, 100,
1053/1062 (1994)
13）R. M. Neal and G. E. Hinton: A view of the EM algorithm that justifies incremental, sparse, and other variants,
Learning in Graphical Models (ed. M. I. Jordan), Kluwer
Academic Publishers, 355/368 (1998)
14）S. Schaal and C. G. Atkeson: Constructive incremental
learning from only local information, Neural Computation,
10, 2047/2084 (1998)15）C. E. Rasmussen, R. M. Neal, G. E. Hinton, D. van Camp,
M. Revow, Z. Ghahramani, R. Kustra and R. Tibshirani:
The DELVE Manual, http://www.cs.toronto.edu/~delve(1996)
16）R. A. Jacobs, M. I. Jordan, S. J. Nowlan and G. E. Hinton:
Adaptive mixtures of local experts, Neural Computation,
3, 79/87 (1991)17）S. R. Waterhouse: Classification and regression using mixtures of experts, Ph. D. Thesis, Department of Engineering,
University of Cambridge (1997)18）M. I. Jordan and R. A. Jacobs: Hierarchical mixtures ofexperts and the EM algorithm, Neural Computation, 6,
181/214 (1994)

本研究は、財団法人テレコム先端技術研究支援センターお《付

よび通信・放送機構の研究委託により実施された．また人工知能研究振興財団の研究支援を受けた．
参 考

文

献

1）L. Xu, M. I. Jordan and G. E. Hinton: An alternativemodel for mixtures of experts, Advances in Neural Information Processing Systems 7 (eds. G. Tesauro, D. S.
Touretzky and T. K. Leen), 633/640, MIT Press (1995)
2）石井 信, 佐藤 雅昭: 正規化ガウス関数ネットワーク，Mixtureof experts と EM アルゴリズム, 日本神経回路学会誌, 6-1,
30/40 (1999)3）S. Ishii and M. Sato: Reconstruction of chaotic dynamics by on-line EM algorithm, Neural Networks, 14-9,
1239/1256 (2001)
4）吉本 潤一郎，石井 信，佐藤 雅昭: 連続力学システムの自動制
御のためのオンライン EM 強化学習法, システム制御情報学会
論文誌, 16-5 (2003)
5）A. P. Dempster, N. M. Laird and D. B. Rubin: Maxi-

A.

確率分布

A. 1

正規分布

録》

p 次元確率ベクトル x が p 次元中心ベクトル µ，p × p 次
逆共分散行列 S をパラメータとする正規分布に従う時，その確率分布 Np (x |µ, S ) は以下で与えられる．

Np (x |µ, S ) ≡ (2π)−p/2 |S|1/2»–
1
× exp − (x − µ)0 S(x − µ)2

(A. 1)(A. 1) 式を x の関数として用いる場合には，これをガウス関数と呼ぶ．

A. 2

Dirichlet 分布
P
確率条件 pi=1 xi = 1, xi ≥ 0 を満足する p 次元確率ベク

計測自動制御学会論文集 第 39 巻 第 5 号

トル x ≡ (x1 , · · · , xp )0 が p 次元ベクトル γ ≡ (γ1 , · · · , γp )02003 年 5 月

511（ b ） VB-M ステップ: Qθ (θ) を以下で更新する．

をパラメータとする Dirichlet 分布に従う時，その確率分布

Dp (x|γ) は以下で与えられる．
Γ(γ1 + · · · + γp + p) γ1γx · · · xpp
Γ(γ1 + 1) · · · Γ(γp + 1) 1
R∞
ここで，Γ(x) ≡ 0 e−t tx−1 dt はガンマ関数である．

Qθ (θ) ≡ Qθ (g)Dp (x|γ) =A. 3

布に従う時，その確率分布 G (x|a, b) は以下で与えられる．

ba a−1x
exp [−bx]Γ(a)Wishart 分布Qθ (µi , Si ) ≡ NN (µi |mi , γi Si )対称行列 ∆ をパラメータとする Wishart 分布に従う時，その確率分布 Wp (X|λ, ∆) は以下で与えられる．
»–
1
Wp (X|λ, ∆) = c · exp − Tr [∆X]
2
ここで，
λ/2

(λ−p−1)/2|∆| |X|Q
2pλ/2 π p(p−1)/4 pn=1 Γ ((λ + 1 − n)/2)c ≡Qθ (Wi , Bi ) ≡DY

NN +1 (wij |vij , βij Ξi )

j=1

×G (βij |γβi /2, γβi λij /2)γ ≡ (γ1 , · · · , γM )0γi ← κT hzi i + γ0ip × p 次正定対称行列 X がスカラー値 λ と p × p 次正定

である．

B.

Qθ (g) ≡ DM (g|γ)×WN (Si |γsi , γsi ∆i )確率変数 x (x ≥ 0) が a と b をパラメータとするガンマ分

A. 4

Qθ (µi , Si )Qθ (Wi , Bi )

i=1

ガンマ分布

G (x|a, b) =M
Y

VB-EMH アルゴリズム

3. 3 節で導出した VB-EMH アルゴリズムは以下のように実装できる．
（ 1 ） 事前分布定数パラメータを設定する．
（ 2 ） ステップ (3a) と (3b) の計算に必要な事後分布 Qθ ,Qξ に関する期待値 h·iθ , h·iξ を適切に初期化する．

γsi ← κT hzi i + γs0
γβi ← κT hzi i + γβ0
κT hzi xi + γ0i m0i
γi
˙
¸
1 “
∆i ←
κT zi xx0 + γ0i m0i m00i
γsi
mi ←

”
+γs0 hσi iξ IN − γi mi m0i(B. 1)˙0¸
Ξi ← κT zi x̃x̃ + hΥi iξ˙¸Vi ← κT zi yx̃0 Ξ−1i“”diag κT hzi yy 0 i + γβ0 hRi iξ − Vi Ξi Vi0Λi ←γβi
ここで，Vi ≡ (vi1 , · · · viD ), Λi ≡ diag (λi1 , · · · , λiD )である．Qθ (θ) に関する期待値 h·iθ は以下で与えられる．

（ 3 ） 自由エネルギー F λ が収束するまで以下のステップ(a)-(c) を繰り返す．
（ a ） VB-E ステップ: QZ を以下で更新する．

QZ (Z) ≡TY

Qz (z(t) = 1)

t=1

exp [Ui (x(t), y(t))]
Qz (zi (t) = 1) ← PM
j=1 exp [Uj (x(t), y(t))]
1 0x hSi iθ x + x0 hSi µi iθ
2¸
1˙1− µ0i Si µi θ − y 0 hBi iθ y
2
2
¸
1 ˙
+y 0 hBi Wi iθ x̃ − x̃0 Wi0 Bi Wi θ x̃2
1
1
+ hlog |Si |iθ + hlog |Bi |iθ22

Ui (x, y) ≡ hlog gi iθ −QZ に関する十分統計量の期待値は以下で求められる．
T

1 X
hzi f (x, u)i =Qz (zi (t) = 1)f (x(t), y(t))
T t=1hlog gi iθ = ψ (γi + 1) − ψMX!
γj + M

j=1

hSi iθ = ∆−1ihSi µi iθ = ∆−1
i mi
˙ 0
¸
µi Si µi θ = m0i ∆−1
i mi + N/γi

(B. 2)

hlog |Si |iθ = − log |∆i | − N log(γsi /2)«„NX
γsi + 1 − n
+
ψ
2
n=1
hBi iθ = Λ−1ihBi Wi iθ = Λ−1
i Vi
˙ 0
¸
−1
Wi Bi Wi θ = Vi0 Λ−1i Vi + DΞihlog |Bi |iθ = − log |Λi | − D log (γβi /2)+Dψ (γβi /2)d
ここで，ψ (x) ≡ dx
Γ(x) はダイガンマ関数である．

512

T. SICE Vol.39 No.5 May 2003

（ c ） VB-H ステップ: Qξ (ξ) を以下で更新する．

Qξ (ξ) ≡MYQξ (σi )Qξ (Υi )Qξ (Ri )

i=1

`
´
−1
Qξ (σi ) ≡ G σi |γσi /2, γσi τσi/2Qξ (Υi ) ≡

N
+1
Y

`
´
−1
G υij |γυi /2, γυi τυij
/2

j=1

Qξ (Ri ) ≡

D
Y

`
´
−1
G ρij |γρi /2, γρi τρij
/2

j=1

γσi ← N γs0 + γσ0 ,γυi ← D + γυ0 ,γρi ← γβ0 + γρ0
`
ˆ
˜
−1
−1 ´
τσi
← γs0 Tr hSi iθ + γσ0 τσ0
/γσi
“`˙
”
¸ ´
−1
0
−1
τυij ←Wi Bi Wi θ j,j + γυ0 τυ0 /γυi
“
”
`
´
−1
−1
τρij
← γβ0 hBi iθ j,j + γρ0 τρ0/γρi
ここで，(·)j,j は行列の j 行 j 列目の要素を表す．Qξ (ξ)
に関する期待値は以下で与えられる．
`´
hΥi iθ = diag τυi1 , · · · , τυi(N +1)

hRi iθ = diag (τρi1 , · · · , τρiD )hσi iθ = τσi

［著 者
吉 本

紹

介］

潤一郎 （正会員）
2002 年 9 月奈良先端科学技術大学院大学情報科学研究科博士後期課程修了．同年 10 月より科学技術振興事業団 CREST 銅谷プロジェクトの研究員となり，現在に至る．博士（工学）．ニュー
ラルネットワーク，統計的学習理論，強化学習の研究に従事．

石 井

信1988 年 3 月東京大学大学院工学系研究科修士課程修了．
（株）リコー中央研究所研究員，ATR
人間情報通信研究所研究員，奈良先端科学技術大学院大学情報科学研究科助教授を経て，現在，同教授．工学博士．非線形力学系，ニューラルネッ
トワーク，強化学習，統計的学習理論，バイオ情報学の研究に従事．

佐 藤

雅

昭
1980 年 3 月大阪大学大学院理学研究科物理学専攻博士課程修了．ニューヨーク大学助手，フロリダ大学助手，ATR 視聴覚機構研究所主任研究員，ATR 人間情報通信研究所主任研究員を経て，現在，ATR 人間情報科学研究所主任研究員．理学博士．非線形力学系，カオス，ニューラルネッ
トワーク，強化学習，ロボット制御，統計的学習理論の研究に従事．