21 離散時間線形モデルの構造特性を活用した改訂RELS 法 常葉大学経営学部紀要 第 2 巻第 2 号，2015 年 2 月，21 － 28 頁
離散時間線形モデルの構造特性を活用した改訂
RELS 法
竹 安 数 博
A Revised Recursive Extended Least Square Method for the Identification of ARMA Model
Kazuhiro TAKEYASU
要 旨
シ ス テ ム 同 定 に お い て、ARMA（Autoregressive Moving Average： 自 己 回 帰 移 動 平 均 ） モ デ ル や、ARMAX （Autoregressive Moving Average with exogenous input：外生入力自己回帰移動平均）モデル等を用いて、過去の入 力信号に基づきシステムパラメータを推定していくやり方が取られる。システムパラメータの不偏推定値を得るアルゴ リズムが各種開発されている。 ここではまず数式モデルを示し、次いで問題と課題を整理し、本論文で改善しようとする点について述べる。
( )
p,
q
次のARMA モデル（自己回帰移動平均モデル）は、p 次の自己回帰過程と q 次の移動平均過程の組み合わせで 表される。 ここで自己回帰過程は、定常エルゴード的正規過程x
( )
t
の標本時系列の線形結合、移動平均過程は、平均値 0、分散 2 eσ
の正規性白色雑音の線形結合である。 なお、ARMA モデルは、一般的には ARMAX モデルにおける制御入力を受けないシステムのモデルと言い換えること ができる。 ARMA モデルにおけるパラメータの推定に際しては、移動平均過程そのものが有色雑音となっており、通常の最小 二乗法を用いてもバイアスのかかった推定値となる。推定値がバイアスを持たないようにするため、拡大最小二乗法、 一般化最小二乗法、逐次拡大最小二乗法、逐次最尤法、補助変数法、Gauss-Newton 法、擬似線形回帰法などがある。 本論文では逐次拡大最小二乗法（Recursive Extended Least Square Method : RELS 法）の改善について述べる。 一般の時系列データは制御入力を受けないシステムである場合も多く、本論文ではその場合についての推定時の計算 時間短縮法等について、次のような視点で改善を図る。RELS 法においてモデルの構造的な前提からくる a priori knowledge を活用することによって、繰り返し各回の計算時間の短縮がなされるとともに収束も早くなることにより、 トータルな計算時間の短縮化が図られる。Abstract
In making system identification, system parameters are estimated using ARMA (Autoregressive Moving Average) model or ARMAX (Autoregressive Moving Average with exogenous input) model based upon the past signal. Many methods which make unbiased estimate have been developed. In ARMA model, AR process is a linear combination of x(t), which is a sample process of a stationary ergodic Gaussian process and MA process is a linear combination of Gaussian noise. Generally, ARMA model is re-stated as the one which does not have control input in ARMAX model. In estimating ARMA model parameters, MA process itself becomes colored noise and we cannot get unbiased estimate by only using Least Square Method. In order to obtain unbiased estimate, such methods as Extended Least Square Method, Generalized Least Square Method, Recursive Extended Least Square Method (RELS Method), Instrumental Variable Method, Gauss-Newton Method, etc. have been developed. In this paper, a revised Recursive Extended Least Square Method is proposed.
In the general time series, there are many cases that they are not affected by the control input, therefore such cases are discussed in this paper. To shorten the calculation time, the following method is proposed. By utilizing the a priori knowledge, which is obtained from the structural presupposition, a new method which shortens the calculation time is proposed. The effectiveness of this method should be examined in various cases.
22 竹 安 数 博
１．はじめに
シ ス テ ム 同 定 に お い て、ARMA（Autoregressive Moving Average： 自 己 回 帰 移 動 平 均 ） モ デ ル や、ARMAX （Autoregressive Moving Average with exogenous input：外生入力自己回帰移動平均）モデル等を用いて、過去の入 力信号に基づきシステムパラメータを推定していくやり方が取られる。システムパラメータの不偏推定値を得るアルゴ リズムが各種開発されている。 ここではまず数式モデルを提示し、次いで問題と課題を整理し、本論文で改善しようとする点について述べる。
( )
p,
q
次のARMA モデル（自己回帰移動平均モデル）は∑
∑
= − = −+
=
+
p i q j j n j n i n i na
x
e
b
e
x
1 1 (1) で表される。 ここで{ }
x
n ：定常エルゴード的正規過程x
( )
t
の標本時系列(
n =1,2,,N,)
{ }
e
n ：平均値 0、分散 2 eσ
の正規性白色雑音 である。( )
( )
q q p pz
b
z
b
z
B
z
a
z
a
z
A
− − − −+
+
+
=
+
+
+
=
1 1 1 11
1
で与えられるA ,( ) ( )
z B z は既約で定常条件、可逆条件、強正実条件を満たすものとする。なお、式(1) は一般的には ARMAX（Autoregressive Moving Average with exogenous input）モデル
q n q n n m n m n p n p n n
e
b
e
b
e
u
c
u
c
x
a
x
a
x
− − − − − −+
+
+
+
+
+
=
+
+
+
1 1 1 1 1 1 (2) における制御入力を受けないシステムのモデルと言い換えることができる。 式(1) のパラメータ{ }
a ,
i{ }
b
j の推定に際しては、右辺そのものが有色雑音となっており、通常の最小二乗法を用いて もバイアスのかかった推定値となる。推定値がバイアスを持たないようにするため、拡大最小二乗法、一般化最小二乗 法、逐次拡大最小二乗法、逐次最尤法、補助変数法、Gauss-Newton 法、擬似線形回帰法などがある [1]-[6]。本論文では逐次拡大最小二乗法（Recursive Extended Least Square Method : RELS 法）の改善について述べる。 一般の時系列データは制御入力を受けないシステムである場合も多く、本論文ではその場合についての推定時の計算 時間短縮法等について、次のような視点で改善を図る。RELS 法においてモデルの構造的な前提からくる a priori knowledge を活用することによって、繰り返し各回の計算時間の短縮がなされるとともに収束も早くなることにより、 トータルな計算時間の短縮化が図られる。計算時間の短縮については、例えば中村・大石[7] では一般化最小二乗法を 用いたものを提案している。そこでは入力、出力およびノイズのある系を対象としている。 上記のような関連するテーマのものはあっても、本論文のようなモデルの構造的な特性に着目したアプローチ・内容 のものは存在しない。以下、２章ではモデルの構造的な特性を明確にするため、３章で述べる準備として、相互相関関 数を活用したARMA モデルのパラメータ推定方法について述べる。３章では RELS 法の改善方法について提案する。 ４章でまとめを示す。
23 離散時間線形モデルの構造特性を活用した改訂RELS 法
２．相互相関関数を活用した ARMAX モデルのパラメータ推定
{ }
x
i の自己相関関数は、定義により[
n n k]
kE
x
x
R
=
+ (3) k kR
R =
− (4) である。また、{ }
x
i と{ }
e
j の相互相関関数は{ }
e
j が白色雑音であるから、l
>
0
に対し q n q n n m n m n p n p n ne
b
e
b
e
u
c
u
c
x
a
x
a
x
− − − − − −+
+
+
+
+
+
=
+
+
+
L
L
L
1 1 1 1 1 1 (2) における制御入力を受けないシステムのモデルと言い換えることができる。 式(1)のパラメータ{
a
i}
,
{ }
b
j の推定に際しては、右辺そのものが有色雑音となっており、通常の最小二 乗法を用いてもバイアスのかかった推定値となる。推定値がバイアスを持たないようにするため、拡大 最小二乗法、一般化最小二乗法、逐次拡大最小二乗法、逐次最尤法、補助変数法、Gauss-Newton 法、擬 似線形回帰法などがある[1]-[6]。本論文では逐次拡大最小二乗法（Recursive Extended Least Square Method : RELS 法）の改善について述 べる。 一般の時系列データは制御入力を受けないシステムである場合も多く、本論文ではその場合について の推定時の計算時間短縮法等について、次のような視点で改善を図る。RELS 法においてモデルの構造的 な前提からくる a priori knowledge を活用することによって、繰り返し各回の計算時間の短縮がなされる とともに収束も早くなることにより、トータルな計算時間の短縮化が図られる。計算時間の短縮につい ては、例えば中村・大石[7]では一般化最小二乗法を用いたものを提案している。そこでは入力、出力お よびノイズのある系を対象としている。 上記のような関連するテーマのものはあっても、本論文のようなモデルの構造的な特性に着目したア プローチ・内容のものは存在しない。以下、２章ではモデルの構造的な特性を明確にするため、３章で 述べる準備として、相互相関関数を活用した ARMA モデルのパラメータ推定方法について述べる。３章 では RELS 法の改善方法について提案する。４章でまとめを示す。
2.
相互相関関数を活用した ARMAX モデルのパラメータ推定
{
i}
x
の自己相関関数は、定義により[
n n k]
kE
x
x
R
+=
(3) k kR
R
=
− (4) である。また、{
i}
x
と{ }
je
の相互相関関数は{ }
je
が白色雑音であるから、l
>
0
に対し( )
[
n n l]
l exl
E
e
x
T
T
=
=
+ (5)( )
=
[
]
=
0
+l n n xel
E
x
e
T
(
−
)
=
[
]
=
0
−l n n exl
E
e
x
T
(
)
[
n n l]
l xel
E
x
e
T
T
−
=
=
− となる。また (5) となる。また[
]
=
l ke
e
E
l
k
l
k
e≠
=
:
0
:
2σ
(6) である。なお、有限個のデータで実際に計算すると、[
n n l]
e
x
E
+(
)
0
>
l
のように本来は無相関で理論的に は 0 であるはずのものが、0 でない値を持つことがある。これを (6) である。なお、有限個のデータで実際に計算すると、E
[
x
ne
n+l]
(
l
>
0
)
のように本来は無相関で理論的には 0 であるは ずのものが、0 でない値を持つことがある。これを( )
[
n n l]
l xel
E
x
e
T
T
− +=
=
(
≅
0
)
(7)(
)
[
n n l]
l exl
E
e
x
T
T
− −=
=
−
(
≅
0
)
[
e
ke
k l]
E
[
e
k le
k]
S
lE
=
=
+ +(
)
0
≅
(8) のように表すことにする。また、上記式(6),(8)の表記に合わせ[
e
k]
eS
oE
=
=
2 2σ
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q
≤
p
となるこ とは周知である[8]。以下、簡単のためp =
q
とする。 式(1)を∑
∑
= = − −+
+
−
=
p i p i n j n j i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
]
T T n T n T p n n p n n nx
x
e
e
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
− − − −L
L
(10)[
]
[
]
T T T T p pb
b
b
a
a
a
b
a
θ
,
,
,
,
,
,
,
,
2 1 2 1=
=
L
L
(11) 式(9)は n n T ne
x
=
θ
Z
+
(12) と表せる。ここで[
]
2 1∑
=−
=
N n n T n Nx
I
θ
Z
(13) を最小にするパラメータ N θˆ を求める。最小二乗推定は∑
∑
= − =⎥
⎦
⎤
⎢
⎣
⎡
=
N n n n N n T n n Nx
1 1 1ˆ
Z
Z
Z
θ
(14) で与えられる。ここでN →∞とすると、式(14)の右辺は各々 ix
と je
との相互相関関数で記述でき、 (7)[
e
ke
k l] [
E
e
k le
k]
S
lE
+=
+=
( )
≅
0
(8) のように表すことにする。また、上記式(6),(8) の表記に合わせ[ ]
e
k eS
oE
2=
σ
2=
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q ≤
p
となることは周知であ る[8]。以下、簡単のためp =
q
とする。 式(1) を∑
∑
= − = −+
+
−
=
p i p i j n j n i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
T]
T n T n T p n n p n n nx
x
e
e
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
−
− −
− (10)( )
[
n n l]
l xel
E
x
e
T
T
− +=
=
(
≅
0
)
(7)(
)
[
n n l]
l exl
E
e
x
T
T
− −=
=
−
(
≅
0
)
[
e
ke
k l]
E
[
e
k le
k]
S
lE
=
=
+ +(
)
0
≅
(8) のように表すことにする。また、上記式(6),(8)の表記に合わせ[
e
k]
eS
oE
=
=
2 2σ
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q
≤
p
となるこ とは周知である[8]。以下、簡単のためp =
q
とする。 式(1)を∑
∑
= = − −+
+
−
=
p i p i n j n j i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
]
T T n T n T p n n p n n ne
e
x
x
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
− − − −L
L
(10)[
]
[
T T]
T T p pb
b
b
a
a
a
b
a
θ
,
,
,
,
,
,
,
,
2 1 2 1=
=
L
L
(11) 式(9)は n n T ne
x
=
θ
Z
+
(12) と表せる。ここで[
]
2 1∑
=−
=
N n n T n Nx
I
θ
Z
(13) を最小にするパラメータ N θˆ を求める。最小二乗推定は∑
∑
= − =⎥
⎦
⎤
⎢
⎣
⎡
=
N n n n N n T n n Nx
1 1 1ˆ
Z
Z
Z
θ
(14) で与えられる。ここでN →∞とすると、式(14)の右辺は各々 ix
と je
との相互相関関数で記述でき、 (11)24 竹 安 数 博 式(9) は
( )
[
n n l]
l xel
E
x
e
T
T
− +=
=
(
≅
0
)
(7)(
)
[
n n l]
l exl
E
e
x
T
T
− −=
=
−
(
≅
0
)
[
e
ke
k l]
E
[
e
k le
k]
S
lE
=
=
+ +(
)
0
≅
(8) のように表すことにする。また、上記式(6),(8)の表記に合わせ[
]
o e kS
e
E
=
=
2 2σ
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q
≤
p
となるこ とは周知である[8]。以下、簡単のためp =
q
とする。 式(1)を∑
∑
= = − −+
+
−
=
p i p i n j n j i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
]
T T n T n T p n n p n n nx
x
e
e
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
− − − −L
L
(10)[
]
[
]
T T T T p pb
b
b
a
a
a
b
a
θ
,
,
,
,
,
,
,
,
2 1 2 1=
=
L
L
(11) 式(9)は n n T ne
x
=
θ
Z
+
(12) と表せる。ここで[
]
2 1∑
=−
=
N n n T n Nx
I
θ
Z
(13) を最小にするパラメータ N θˆ を求める。最小二乗推定は∑
∑
= − =⎥
⎦
⎤
⎢
⎣
⎡
=
N n n n N n T n n Nx
1 1 1ˆ
Z
Z
Z
θ
(14) で与えられる。ここでN→∞とすると、式(14)の右辺は各々 ix
と je
との相互相関関数で記述でき、 (12) と表せる。ここで( )
[
n n l]
l xel
E
x
e
T
T
− +=
=
(
≅
0
)
(7)(
)
[
n n l]
l exl
E
e
x
T
T
− −=
=
−
(
≅
0
)
[
e
ke
k l]
E
[
e
k le
k]
S
lE
=
=
+ +(
)
0
≅
(8) のように表すことにする。また、上記式(6),(8)の表記に合わせ[
]
o e kS
e
E
=
=
2 2σ
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q
≤
p
となるこ とは周知である[8]。以下、簡単のためp =
q
とする。 式(1)を∑
∑
= = − −+
+
−
=
p i p i n j n j i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
]
T T n T n T p n n p n n nx
x
e
e
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
− − − −L
L
(10)[
]
[
]
T T T T p pb
b
b
a
a
a
b
a
θ
,
,
,
,
,
,
,
,
2 1 2 1=
=
L
L
(11) 式(9)は n n T ne
x
=
θ
Z
+
(12) と表せる。ここで[
]
2 1∑
=−
=
N n n T n Nx
I
θ
Z
(13) を最小にするパラメータ N θˆ を求める。最小二乗推定は∑
∑
= − =⎥
⎦
⎤
⎢
⎣
⎡
=
N n n n N n T n n Nx
1 1 1ˆ
Z
Z
Z
θ
(14) で与えられる。ここでN→∞とすると、式(14)の右辺は各々 ix
と je
との相互相関関数で記述でき、 (13) を最小にするパラメータ( )
[
]
l l n n xeT
e
x
E
l
T
− +=
=
(
≅
0
)
(7)(
)
[
]
l l n n exT
x
e
E
l
T
− −=
=
−
(
≅
0
)
[
]
[
]
l k l k l k kS
e
e
E
e
e
E
=
=
+ +(
)
0
≅
(8) のように表すことにする。また、上記式(6),(8)の表記に合わせ[
e
k]
eS
oE
=
=
2 2σ
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q
≤
p
となるこ とは周知である[8]。以下、簡単のためp =
q
とする。 式(1)を∑
∑
= = − −+
+
−
=
p i p i n j n j i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
T]
T n T n T p n n p n n nx
x
e
e
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
− − − −L
L
(10)[
]
[
]
T T T T p pb
b
b
a
a
a
b
a
θ
,
,
,
,
,
,
,
,
2 1 2 1=
=
L
L
(11) 式(9)は n n T ne
x
=
θ
Z
+
(12) と表せる。ここで[
]
2 1∑
=−
=
N n n T n Nx
I
θ
Z
(13) を最小にするパラメータ N θˆ を求める。最小二乗推定は∑
∑
= − =⎥
⎦
⎤
⎢
⎣
⎡
=
N n n n N n T n n Nx
1 1 1ˆ
Z
Z
Z
θ
(14) で与えられる。ここでN→∞とすると、式(14)の右辺は各々 ix
と je
との相互相関関数で記述でき、 を求める。最小二乗推定は( )
[
n n l]
l xel
E
x
e
T
T
− +=
=
(
≅
0
)
(7)(
)
[
n n l]
l exl
E
e
x
T
T
− −=
=
−
(
≅
0
)
[
e
ke
k l]
E
[
e
k le
k]
S
lE
=
=
+ +(
)
0
≅
(8) のように表すことにする。また、上記式(6),(8)の表記に合わせ[
e
k]
eS
oE
=
=
2 2σ
と表すことにする。 現代制御理論において、系の状態空間表現を行い、ARMA 過程表現に落とし込むと、q
≤
p
となるこ とは周知である[8]。以下、簡単のためp =
q
とする。 式(1)を∑
∑
= = − −+
+
−
=
p i p i n j n j i n i na
x
b
e
e
x
1 1 (9) と書き直し、ベクトルを以下のように定義する。[
]
[
]
T T n T n T p n n p n n nx
x
e
e
e
x
Z
,
,
,
,
,
,
1 1−
=
−
−
=
− − − −L
L
(10)[
]
[
]
T T T T p pb
b
b
a
a
a
b
a
θ
,
,
,
,
,
,
,
,
2 1 2 1=
=
L
L
(11) 式(9)は n n T ne
x
=
θ
Z
+
(12) と表せる。ここで[
]
2 1∑
=−
=
N n n T n Nx
I
θ
Z
(13) を最小にするパラメータ N θˆ を求める。最小二乗推定は∑
∑
= − =⎥
⎦
⎤
⎢
⎣
⎡
=
N n n n N n T n n Nx
1 1 1ˆ
Z
Z
Z
θ
(14) で与えられる。ここでN→∞とすると、式(14)の右辺は各々 ix
と je
との相互相関関数で記述でき、 (14) で与えられる。ここでN →∞とすると、式(14) の右辺は各々x
iとe
jとの相互相関関数で記述でき、
−
−
−
=
t
r
I
T
T
R
b
a
2 e Tσ
(15) となる。ここで
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
R
R
R
R
R
R
R
R
R
p p p p
R
=
− − 0 2 0 1 1 0,
,
0
,
0
,
,
,
0
,
,
,
T
T
T
T
T
T
p p
T
[
]
T pR
R
R
1,
2,
,
=
r
[
]
T pT
T
T
1,
2,
,
=
t
である。式(15) を書き直すと、⎥
⎦
⎤
⎢
⎣
⎡−
⎥
⎦
⎤
⎢
⎣
⎡
−
−
=
⎥
⎦
⎤
⎢
⎣
⎡
t
r
I
T
T
R
b
a
2 e Tσ
(15) となる。ここで⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
R
R
R
R
R
R
R
R
R
p p p pL
M
O
M
M
L
L
R
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
− − 0 2 0 1 1 0,
,
0
,
0
,
,
,
0
,
,
,
T
T
T
T
T
T
p pL
M
O
M
M
L
L
T
[
R
,
R
,
,
R
p]
T 2 1L
=
r
[
T
,
T
,
,
T
p]
T 2 1L
=
t
である。式(15)を書き直すと、r
Tb
Ra
−
=
−
(16)t
b
a
T
+
=
−
2 e Tσ
(17) となる。式(17)より、(
t
T
a
)
b
T e+
=
21
σ
(18) が得られるが、これを式(16)へ代入することにより⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
=
−r
Tt
TT
R
a
2 1 21
1
e T eσ
σ
(19) が得られる。これを式(18)に代入してb
が求められる。パラメータ推定の全体のアルゴリズムは、表 1 のようなものが考えられ、収束の判定は ( ) ( ) ( ) ( ) ( ) ( ) ( )∑
= − − − −⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
+
−
=
p i l i l i l i l i l i l i lb
b
b
a
a
a
D
1 1 1 1 1ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
(20) 表 1. パラメータ推定の全体アルゴリズム (16)t
b
a
T
+
=
−
2 e Tσ
(17) となる。式(17) より、(
t
T
a
)
b
T e+
=
1
2σ
(18) が得られるが、これを式(16) へ代入することにより⎥
⎦
⎤
⎢
⎣
⎡−
⎥
⎦
⎤
⎢
⎣
⎡
−
−
=
⎥
⎦
⎤
⎢
⎣
⎡
t
r
I
T
T
R
b
a
2 e Tσ
(15) となる。ここで⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
R
R
R
R
R
R
R
R
R
p p p pL
M
O
M
M
L
L
R
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
− − 0 2 0 1 1 0,
,
0
,
0
,
,
,
0
,
,
,
T
T
T
T
T
T
p pL
M
O
M
M
L
L
T
[
]
T pR
R
R
,
,
,
2 1L
=
r
[
]
T pT
T
T
,
,
,
2 1L
=
t
である。式(15)を書き直すと、r
Tb
Ra
−
=
−
(16)t
b
a
T
=
+
−
2 e Tσ
(17) となる。式(17)より、(
t
T
a
)
b
T e+
=
21
σ
(18) が得られるが、これを式(16)へ代入することにより⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
=
−r
Tt
TT
R
a
2 1 21
1
e T eσ
σ
(19) が得られる。これを式(18)に代入してb
が求められる。パラメータ推定の全体のアルゴリズムは、表 1 のようなものが考えられ、収束の判定は ( ) ( ) ( ) ( ) ( ) ( ) ( )∑
= − − − −⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
+
−
=
p i l i l i l i l i l i l i lb
b
b
a
a
a
D
1 1 1 1 1ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
(20) 表 1. パラメータ推定の全体アルゴリズム (19) が得られる。これを式(18) に代入してb
が求められる。パラメータ推定の全体のアルゴリズムは、表 1 のようなもの が考えられ、収束の判定は25 離散時間線形モデルの構造特性を活用した改訂RELS 法 ( ) ( ) ( ) ( ) ( ) ( ) ( )
∑
= − − − −
−
+
−
=
p i il l i l i l i l i l i lb
b
b
a
a
a
D
1 1 1 1 1ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
(20) 表 1. パラメータ推定の全体アルゴリズム ステップ 1 : 観測データ{ }
x
n から自己相関関数{ }
Rˆ
k を計算する ステップ 2 : 正規性白色雑音を発生させ、{ }
e
n の初期値とする ステップ 3 :{ }
Tˆ
k を計算する ステップ 4 : 式 (18),(19) よりa
ˆ
( ) ( )l,
b
ˆ
l を推定する ステップ 5 : 式 (9) より{ }
e
n を推定する ステップ 7 : ステップ 3,4,5 をa
ˆ
( ) ( )l,
b
ˆ
l が収束するまで繰り返す において ( )l<
ε
D
(21) がm
回連続で成立したときとする。T
は上三角行列のため、その性質を利用して計算時間の短縮が可能となる。従来の方法であれば、T
の全成分が埋まっ た形でパラメータ推定計算がなされるが、ここでは構造的な前提からくるa priori knowledge を活用することによっ て ステップ 1 : 観測データ{
n}
x
から自己相関関数{
}
kR
ˆ
を計算する ステップ 2 : 正規性白色雑音を発生させ、{
}
ne
の初期値とする ステップ 3 :{ }
T
kˆ
を計算する ステップ 4 : 式(18),(19)より ( )l ( )lb
a
ˆ
,
ˆ
を推定する ステップ 5 : 式(9)より{
}
ne
を推定する ステップ 7 : ステップ 3,4,5 を ( )l ( )lb
a
ˆ
,
ˆ
が収束するまで繰り返す において ( )ε
<
lD
(21) がm
回連続で成立したときとする。T
は上三角行列のため、その性質を利用して計算時間の短縮が可能となる。従来の方法であれば、T
の全成分が埋まった形でパラメータ推定計算がなされるが、ここでは構造的な前提からくる a priori knowledgeを活用することによって 解のロバスト性 計算時間の短縮 が期待できる。なお、式(7),(8)を用いた場合、式(15)に相当するものは次のようになる。⎥
⎦
⎤
⎢
⎣
⎡−
⎥
⎦
⎤
⎢
⎣
⎡
−
−
=
⎥
⎦
⎤
⎢
⎣
⎡
−t
r
S
T
T
R
b
a
1 T (22) となる。ただし、⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
R
R
R
R
R
R
R
R
R
p p p pL
M
O
M
M
L
L
R
( ) ( )⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
T
T
T
T
T
T
T
T
T
p p p pL
M
O
M
M
L
L
T
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
S
S
S
S
S
S
S
S
S
p p p pL
M
O
M
M
L
L
S
解のロバスト性 ステップ 1 : 観測データ{
}
nx
から自己相関関数{
}
kR
ˆ
を計算する ステップ 2 : 正規性白色雑音を発生させ、{
n}
e
の初期値とする ステップ 3 :{ }
T
kˆ
を計算する ステップ 4 : 式(18),(19)より ( )l ( )lb
a
ˆ
,
ˆ
を推定する ステップ 5 : 式(9)より{
n}
e
を推定する ステップ 7 : ステップ 3,4,5 を ( )l ( )lb
a
ˆ
,
ˆ
が収束するまで繰り返す において ( )ε
<
lD
(21) がm
回連続で成立したときとする。T
は上三角行列のため、その性質を利用して計算時間の短縮が可能となる。従来の方法であれば、T
の全成分が埋まった形でパラメータ推定計算がなされるが、ここでは構造的な前提からくる a priori knowledgeを活用することによって 解のロバスト性 計算時間の短縮 が期待できる。なお、式(7),(8)を用いた場合、式(15)に相当するものは次のようになる。⎥
⎦
⎤
⎢
⎣
⎡−
⎥
⎦
⎤
⎢
⎣
⎡
−
−
=
⎥
⎦
⎤
⎢
⎣
⎡
−t
r
S
T
T
R
b
a
1 T (22) となる。ただし、⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
R
R
R
R
R
R
R
R
R
p p p pL
M
O
M
M
L
L
R
( ) ( )⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
T
T
T
T
T
T
T
T
T
p p p pL
M
O
M
M
L
L
T
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
S
S
S
S
S
S
S
S
S
p p p pL
M
O
M
M
L
L
S
計算時間の短縮 が期待できる。なお、式(7),(8) を用いた場合、式 (15) に相当するものは次のようになる。
−
−
−
=
−t
r
S
T
T
R
b
a
1 T (22) となる。ただし、
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
R
R
R
R
R
R
R
R
R
p p p p
R
( ) ( )
=
− − − − − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
T
T
T
T
T
T
T
T
T
p p p p
T
=
− − − − 0 2 1 2 0 1 1 1 0,
,
,
,
,
,
,
,
,
S
S
S
S
S
S
S
S
S
p p p p
S
26 竹 安 数 博
