関連ベクターマシンに基づく非線形回帰モデルと変化点探索
九州大学大学院数理学府
九州大学大学院数理学研究院

1立石 正平小西 貞則

はじめに回帰分析において，複雑な非線形構造を内包する現象を解明するために，スプライン関数やカーネル関数に基づく非線形回帰モデルが用いられる．しかし，ある点を境に平均構造が突然変化するデータに対しては，非線形回帰モデルの当てはめが困難となり，データから有効な情報を得ることが難しくなる．このようなデータに対しては，変化点を適切に捉えることが重要となる (井元・小西 (2000))．
本報告では，Tipping (2001) による関連ベクターマシン (RVM: Relevance Vector Machine) を適用することで，変化点をもつ非線形構造を適切に推定する方法を提案する．基
底関数展開に基づく非線形回帰モデルを仮定する場合，調整パラメータの選択が本質的となる．この問題に対して，モデル評価基準 (Konishi and Kitagawa, 2008) を導出し，調整パラメータの選択を行う．最後に，提案手法の有効性を数値例を通して検証する．

2

基底関数展開に基づく非線形回帰モデル目的変数 Y と説明変数 X に関して，n 組のデータ {(yi , xi ); i = 1, 2, · · · , n} が観測されたとする．このとき，説明変数と目的変数の間に次の関係が成り立つと仮定する．

yi = u(xi ; w) + i ,(i = 1, 2, · · · , n) .

(1)

ここで， は誤差項，u(x; w) は，未知のパラメータベクトル w に依存する回帰関数である．
次に，以下のように回帰関数として基底関数と呼ばれる既知の非線形関数 b(x) に基づくモデルを考える．

u(x; w) = w0 b(x) .

(2)

ただし，b(x)=(1,b1 (x),b2 (x),· · · ,bm (x))0 は (m + 1) 次元基底関数ベクトル，w = (w0 , w1 ,· · · , wm )0 は (m + 1) 次元係数パラメータベクトルとする．誤差項 i は互いに独立に平均0，分散 σ 2 の正規分布に従うとすると，基底関数に基づく回帰モデルは次のように表される．
{}(yi − w0 b(xi ))222 − 21.
(3)f (yi |w, σ ) = (2πσ ) exp −2σ 21

3

変化点探索と曲線推定

本報告では，変化点の探索と曲線推定の 2 段階の推定を行う．まず，第 1 段階では，変化点を捉えるための基底関数に基づく非線形回帰モデルを仮定する．次に，第 2 段階では，第 1 段階で用いた基底関数に加え，さらに，ガウス型基底関数を融合した新たな非線形回帰モデルを仮定する．

3.1

関連ベクターマシンに基づく変化点探索

n 組のデータ {(yi , xi ); i = 1, 2, · · · , n} が観測されたとき，観測値ベクトル y = (y1 , · · · , yn )0の分布は (3) 式より次の分布に従う．
y|w, β ∼ N (Φw, σ 2 In )(4)

3.1 節では，β −1 = σ 2 とし，β はハイパーパラメータとして扱う．Φ は基底関数からなるn × (n + 1) 行列で，
1 φ1 (x1 ) · · · φn (x1 ) 1 φ1 (x2 ) · · · φn (x2 ) 

Φ= .
(5)

 ..

1 φ1 (xn ) · · · φn (xn )とする．ここで，データ数 n と同じ数の基底関数を設定することが重要となる．変化点探索における基底関数 φ(x) として，次の関数を用いる．
()
||x − xj ||2exp, (x < xj )h2
φj (x) =(6)0, () (x = xj ) , j = 1, 2, · · · , n.
2
||x − xj || − exp, (x > xj )h2
図 1 には基底関数 φ(x) の例を与える．φ(x) は，その中心の点の前後で値が大きく変化する関数であるので，係数パラメータの推定値の絶対値が大きくなる点が変化点である可能性があると考えられる．
次に係数パラメータ w に次の事前分布を設定する．

w|α ∼ N (0, A−1 ).

(7)

ただし，α = (α0 , · · · , αn )0 ，A = diag(α0 , · · · , αn ) である．このような事前分布は，ARD(関連自動決定) 事前分布と呼ばれる．このとき，ベイズの定理より，w の事後分布は

w|y, α, β ∼ N (µ, Σ)(8)Σ = (βΦ0 Φ + A)−1 ,(9)µ = βΣΦ0 y(10)となる．ただし，

21.00.50.0−0.5−1.00.00.20.40.60.81.0x

図 1: 基底関数 φ(x) の構成例．とする．係数パラメータベクトル w の推定は，(10) 式の事後平均ベクトルで行う．
ここで，ハイパーパラメータ α，β の推定について考える．観測値ベクトル y の周辺尤度関数は∫p(y|α, β) = f (y|w, β)π(w|α)dw(11)であり，自然共役の関係性から，

y|α, β ∼ N (0, β −1 In + ΦA−1 Φ0 )(12)となる．この周辺尤度を最大にするハイパーパラメータの推定量 α̂，̂β は次で与えられる．

γj||y − Φµ||2−1α̂j = 2 , β̂ =∑µjn − n+1j=0 γj(13)ただし，γj = 1 − αj Σjj である．α̂j と β̂ は互いに依存しているため，上式で繰り返し更新を行う．
以上のように，周辺尤度最大化に基づいて係数パラメータの推定とハイパーパラメータの推定を逐次計算によって行う手法は関連ベクターマシン (RVM: Relevance VectorMachine) と呼ばれ，Tipping (2001) によって提案された．RVM は，推定に不要な基底関数にかかる係数パラメータを真に 0 に推定する．従って，切片項を除く n 個の係数パラメータの中で，0 と推定されなかった基底関数に対応する点が変化点の候補と考えることができ，変化点の可能性がある点を {τ1 , τ2 , · · · , τt } と絞り込むことができる．ここで，変化点候補の組み合わせの集合族 {∅, τ1 , · · · , (τ1 , τ2 ), · · · , (τ1 , τ2 , · · · , τt )} の元の 1 つを T = {t1 , · · · , tnT } (nT < n) とする．

3

3.2

正則化法に基づく曲線推定

3.1 節では，(6) 式の基底関数 φ(x) を用いた非線形回帰モデルを仮定し，RVM を適用することで，変化点探索を行った．次に，(3) 式の基底関数ベクトル b(x) をb(x) = (1, ψ1 (x), · · · , ψm (x), φ̃T (x)0 )0(14)とした基底関数展開に基づく非線形回帰モデルを仮定する．ここで，φ̃T (x) = (φt1 (x), · · ·, φtnT (x))0 とする．ψ(x) は次で与えられるガウス型基底関数である．
)(||x − ck ||2, k = 1, 2, · · · , m.
(15)
ψk (x) = exp2νs2k
ただし，m (< n) はガウス型基底関数の個数，ck は関数の中心，s2k は関数の広がりを表す量である．また，ν (> 0) は ψk (x) の広がりの程度を調整するパラメータである (Ando etal., 2008)．
ここで，ck と s2k の決定法を述べる．n 個のデータ {x1 , · · · , xn } をクラスタリング手法の一つである K-means 法によって m 個のクラスタ C1 , · · · , Cm に分割し，各クラスタ Ckに含まれる nk 個のデータに基づいて ck と s2k は，

ck =1 ∑1 ∑xi , s2k =||xi − ck ||2nk x ∈Cnk x ∈Ciik(16)kで決定する．
複雑な現象を近似するために非線形モデルを当てはめるとき，モデルのパラメータを最尤法によって推定すると，しばしばデータに過度に依存したモデルが推定される．そこ
で，u(x; w) の滑らかさが失われるにつれて増加する罰則項 (正則化項) を対数尤度関数に付け加えたn∑nλ2lλ (w, σ ) =log f (yα |w, σ 2 ) −H(w)2α=1の最大化，すなわち正則化法によって推定する．ここで，λ (> 0) は平滑化パラメータと呼ばれ，モデルの当てはまりのよさと曲線の滑らかさを制御するハイパーパラメータである．また，H(w) は正則化項と呼ばれ，w の 2 次形式 w 0 Kw などが用いられる．ただし
K は既知の m + nT + 1 次の非負定値行列である．このモデルのパラメータ wとσ 2 の正則化最尤推定量は次で与えられる．

ŵ = (B 0 B + nλσ̂ 2 K)−1 B 0 y,σ̂ 2 =1(y − B ŵ)0 (y − B ŵ) .
n

ただし，

 
1 ψ1 (x1 ) ψ2 (x1 )b(x1 )0 b(x2 )0   1 ψ1 (x2 ) ψ2 (x2 ) 
B= =  ..
..
..
.
..
  .
.
.
0
1 ψ1 (xn ) ψ2 (xn )b(xn )4ψm (x1 ) φ̃T (x1 )0ψm (x2 ) φ̃T (x2 )0 


..
..
.
.
· · · ψm (xn ) φ̃T (xn )0······...

である．これらの推定量は，分散 σ 2 の推定を通して互いに依存しているので，以下のようなアルゴリズムで計算できる．

[step1] 初期値を σ 2(0) とし，λ をあらかじめ与えておく．k = 1 とする．
(k)0[step2] w = (B B + nλσ 2(k−1) K)−1 B 0 y で更新する．
[step3] σ 2(k) = n1 (y − Bw(k) )0 (y − Bw(k) ) を求め，十分小さな δ に対して，(σ 2(k) −σ 2(k−1) )2 < δ を評価する．この条件を満たすまで step2 を繰り返し，条件を満たしたパラメータを推定値とする．

4

モデル評価基準

前節で述べた推定したモデルは，ガウス型基底関数の個数 m，平滑化パラメータ λ，ガウス基底関数の幅調整パラメータ ν に依存しており，この選択には以下のモデル評価基準を用いる．
一般化情報量規準 GIC (Konishi and Kitagawa, 1996)
データを発生した真の分布 G(x) から得られた n 個 のデータ x = {x1 , x2 , · · · , xn } に基づく経験分布関数を Ĝ(x) とするとき，統計的汎関数 T (·) によって定義される一般の推定量 θ̂ = T (Ĝ) = (T1 (Ĝ), T2 (Ĝ), · · · , Tp (Ĝ))0 に基づくモデル f (x|θ̂) の評価基準として，次で与えられる一般化情報量規準 GIC (Konishi and Kitagawa (1996)) がある．
{}
nn∑2∑∂logf(x|θ)
αGIC = −2log f (xα |θ̂) +tr T (1) (xα ; Ĝ).
n∂θ 0ˆ
α=1α=1θ
ただし，T (1) (xα ; Ĝ) = (T1 (xα ; Ĝ), · · · , Tp (xα ; Ĝ))0 は経験影響関数とよばれ，その第 i(1)
要素 Ti (xα ; Ĝ) は次で定義される．
(1)

(1)

Ti ((1 − )Ĝ + δα ) − Ti (Ĝ)
.
→0

(1)

Ti (xα ; Ĝ) = lim

ここで，δα は点 xα 上に確率 1 をもつ分布関数である．この考えを，前述の正則化最尤法によって推定された基底関数展開に基づく非線形モデルに適用すると，回帰モデルにおける GIC は

GIC = n(log 2π + 1) + n log σ̂ 2 + 2tr{R−1 Q}(17)で与えられる．ただし，

1 01  B B + nλσ̂ K σ̂ 2 B Λ1n R = ,n1 0nσ̂ 2ΛB1σ̂ 2 n2σ̂ 2
1 0 21 0 31 001  σ̂ 2 B Λ B − λK ŵ1n ΛB 2σ̂ 4 B Λ 1n − 2σ̂ 2 B Λ1n 
Q =
 1

1 01 0 4n0 3nσ̂ 2Λ
B−ΛB
Λ
1−111nnnn2σ̂ 42σ̂ 24σ̂ 64σ̂ 2025(18)(19)

であり，1n は要素が全て 1 の n 次元ベクトル，Λ は

Λ = diag(y1 − ŵ0 b(x1 ), · · · , yn − ŵ0 b(xn ))とする．

5適用例

データはガウスノイズに対して N (0, η 2 ) を仮定した次の式に従って，100 個のデータを発生させた．

yi = sin (2πxi ) + I(xi ≥ x50 ) + i ,(i = 1, · · · , n)(20)

ただし，I(xi ≥ a) = 0 (xi < a)，I(xi ≥ a) = 1 (xi ≥ a) である．ここで，データ x は[0, 1] を 100 等分した点とし，ノイズの標準偏差 η は，η = 0.2 とした．変化点は 50 番目のデータ点 x50 = 0.495 に設定した．
(20) 式を真の曲線とするデータに対して提案手法を適用した結果を図 2 に示す．今回
は，(6) 式の基底関数の幅は十分大きい値で固定し，RVM で 0 と推定されなかった係数の中で絶対値の値が上位 2 つの係数と対応する点を変化点の候補とした．変化点の位置，ガウス型基底関数の個数 m，平滑化パラメータ λ，ガウス基底関数の幅調整パラメータ ν は一般化情報量規準 GIC で選択した．
100 回の実験を行った結果，変化点の数を正確に選択できた回数は 88 回であり，その際の推定された変化点の平均値は 0.491，標準偏差は 7.74 × 10−3 であった．

5.1

数値実験

5.2

ナイルデータ1871 年から 1970 年までのナイル河の年間河川流量 (108 m3 ) を計測したものである．
cobb(1978) は，1898 年を境に年間河川流量が突然変化していることを示している．
このナイルデータに対して，提案手法を適用した結果を図 3 に示す．(6) 式の基底関数の幅は十分大きい値で固定し，RVM で 0 と推定されなかった係数の中で絶対値の値が上位 2 つの係数と対応する点を変化点の候補とした．変化点の位置，ガウス型基底関数の個数 m，平滑化パラメータ λ，ガウス基底関数の幅調整パラメータ ν は一般化情報量規準
GIC で選択した結果，1898 年が変化点として選択された．

5.3

硬貨データ

. 0.0254mm) を計測し1945 年から 1989 年までのリンカーン 1 セント硬貨の厚さ (1mil =.
たものである (Scott, 1992)．硬貨の厚さは第 2 次世界大戦中に減らされ，1960 年ごろに元6

1.00.00.5

y1.00.00.5

y0.00.20.40.60.8

1.00.00.20.4x0.60.8

1.0x5253

600

54800

55yy1000

5657

1200581400

図 2: 人工データ解析例．左図は，実線が真の曲線．右図は，実線が推定曲線を示す．

18801900

19201940

19501960

1960

1970

1980

1990

xx

図 3: ナイルデータ．実線が推定曲線を示す．図 4: 硬貨データ．実線が推定曲線を示す．の厚さに戻されたが，1970 年代に再び減らされたことを示している．Speckman(1994) やGijbels and Goderniaux(2004) では，1958 年と 1974 年に変化が起こったと推定している．
この硬貨データに対して，提案手法を適用した結果を図 4 に示す．(6) 式の基底関数の幅は十分大きい値で固定し，RVM で 0 と推定されなかった係数の中で絶対値の値が上位3 つの係数と対応する点を変化点の候補とした．変化点の位置，ガウス型基底関数の個数m，平滑化パラメータ λ，ガウス基底関数の幅調整パラメータ ν は一般化情報量規準 GICで選択した結果，1958 年と 1974 年が変化点として選択された．

76

まとめ

基底関数展開に基づく非線形回帰モデルを仮定し，関連ベクターマシンを適用することで，変化点をもつ非線形構造を適切に推定する方法を提案した．また，変化点の位置，ガウス型基底関数の個数と調整パラメータ，そして平滑化パラメータの選択に情報量規準
GIC を用いた．そして，平均構造に変化点を含む実データと人工データを用いた数値例を通して，提案手法の有効性を示した．

参考文献[1] 井元清哉，小西貞則 (2000). B-スプラインによる不連続曲線推定と変化点探索. 研究集会”多変量解析と漸近理論に関する新展開”. 広島大学．
[2] Ando, T., Konishi, S. and Imoto, S. (2008). Nonlinear regression modeling via regularized radial basis function networks. Journal of Statistical Planning and Inference.
138, 3616–3633.
[3] Cobb, G. W. (1978). The problem of the Nile: Conditional solution to a changepoint
problem. Biometrika. 65, 243–251.
[4] Gijbels, I. and Goderniaux, A.-C. (2004). Bandwidth selection for change point estimation in nonparametric regression. Technometrics, 46, 76–86.
[5] Konishi, S. and Kitagawa, G. (1996). Generalised information criteria in model selection. Biometrika. 83, 875–890.
[6] Konishi, S. and Kitagawa, G. (2008). Information Criteria and Statistical Modeling.
Springer.
[7] Scott, D. W. (1992). Multivariate density estimation. Theory, practice and visualization. New York: Wiley.
[8] Speckman, P. L. (1997). Detection of change-points in nonparametric regression.
Unpublished manuscript.
[9] Tipping, M.E. (2001). Sparse Bayesian learning and the relevance vector machine.
Journal of Machine Learning Research. 1, 211–244.

8