多次元ARモデルによるシステム解析

多次元 AR モデノレによるシステム解析

石黒

真木夫

''''""""""""111""""'''""111"111"""111"""111'111""""111"'"""111"""'''""""111"111"111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111

はじめに

システム解析家の道具箱に入れておいてもらいたい方 法を提案しようと思う.多次元 AR モテ'ルを「分析する」 ための道具である. 身もふたもない言い方をすれば， r赤池と中 111 (1972) によって提案された“時系列の解析と制御のためのプロ グラムパッケージ TIMSAC" に含まれる MULNOS の機能を拡張して応用の範囲を広げるとともに，対話型 にして使いがってを良くしたフ。ログラム J の提案という ことになる.将来の参照の際の便利のためにこれを AR. DOCK と名つeけておく. rAR モデルドック j のつもり である.意味は「人間ドック」から類推されたい. ARDOCK の特徴は，以下のとおりである. -システムのコンポーネント聞の信号経路を切断するこ との影響をモデル上で調べることができる.結果を見な がら対話的に切断，接続を試みることができる. ・ノイズ源相互が完全には相互に無相関でない場合にも 対処できる. 2 番目の主張は， MULNOS を使ったことのない読者 には良くわからないかもしれないが， MULNOS を利用 するに当つての 1 つの障害がある程度回避できるという ことである.

2

.

多次元 AR 毛デルとノイズ寄与率

時系列モデルの使われる分野は非常に広く，工学的な データや，自然科学の分野のデータはもちろんのこと， いわゆる文科系のデータ解析においても効果をあげてい る.その中でとりわけ愛用されているのが AR モデルで ある. 2.1 モデルのあてはめ 多次元 AR モデルは いしぐろ まきお文部省統計数理研究所 干 106 港区南麻布 4-6-7 1989 年 10 月号

x (i)= 主 A(m)x(i -m)+õ(i)

(1)

η1. =1 で表わされる.ここで， {x(i)} は h 次元ベクトルの系列， {A(m)} は h 次元係数行列である.時系列 j{x(i)} の振舞 いは，次数M， 係数行列 {A(m) : i=I ， 2 ， … ， M} および ノイズ系列 {õ(i)} の分散共分散行列 Z などのパラメータ によって決定される. 与えられたデータ {x (i)

:

i=I ， 2 ， … ， N} が(1)式の機構 で生成されると想定すると，最尤法でパラメータを推定 することができる.データが十分に長い場合， {A(m): i=I ， 2 ， … ， M} の推定値が満たすべき最尤方程式は近似 的に N M N

L

:

x(i)xT(i-m)=

L

:

A(k)

L

:

x

(

i

-k)xT(i-m

i=1 .¥:=1 i=1 m=I,

2

,"',

M

ω) となる.ここで (T) は転置の意味である. Z に関してなにも制約を置かない場合には， の推定値が定まったという条件のもとで {A(m)}

(

i

)

=x(i) ーさ A(m)x(i-m)(i=1 ， 2， … ， N) ，

N 拙 =1 B =

L

:

(i)T(i) として，

L:=長B

が最尤法による Z の推定値になる. Z としてブロック対角型の構造，たとえば

L

:

=

(

L

:

pp

_...,

¥

L

:

qq/ を仮定すると，

B ー (l!.pp

l

!

?

q -¥Bqp Bqq) と B を分割して得られる {Bpq} を使って

L:pp =占Bpp， L:qq =会Bqq

(3) (4) と推定されることになる.対角プロックが一般に L 個の 場合にも同様である. Z としてプロック対角裂の構造を仮定するモデルの A

1

C'主， プロックの個数を L として (35)

5

4

7

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

AIC= 去 Nkt

l

o

g

27r+ N

l

o

g

det

L

:

!!+Nkt +(k_t十 l)kt (5) で与えられる.

2

.

2

パワースペクトルとノイズ寄与率 多次元 AR モデルのパラメータが推定できると，周波 数応答関数

A(f)=[I-EA(m)rwm]→ (O~f云0.5)

(6) 間 =1 と，クロスパワースベクトル， P(f)=A(f)

L

:

A*

(

f

)

(O~f云 0.5)η がただちに計算できる.ここで I は単位行列. (*)は共 役転置の意味である. Z が(4)式のようなプロック構造をもっていると. (竹式 から (Pvv

(

f

)

Ppq

(

f

)¥

(App(

f)

A阿(f)\ P(f)

=

(

:PP

,':.:

-DPq;J..:). A(f)

=

(

-;PP,'-:,:

7

q;J..: ) ¥Pqp(f) Pqq(

f)

J ¥Aqp

(

f

)

Aqq(f)J として，

P pp(f) =App

(

f

)

L

:

ppApp*(f) +Apq(f)

L

:

qqApq*(f) Pqq(f) =Aqp(f)

L

:

ppAqp*(f) +Aqq(f) L: qqAqq本 (f)

(

8

)

が得られる. ベクトル x (i) の第 j 成分のパワ}スベクトルは行列 P(f) の(j.j) 成分として推定されるが. (8)式の意味する ところは . Pjj が j の値によって P jj(f)= [App (f)L: ppApp本 (f)Jjj + [Apq(f)

L

:

qqApq*

(f

)Jjj あるいは P jj(f) =[Aqp(f)

L

:

ppAqp*(

f

)

Jjj + [Aqq(f)

L

:

qqAqq*(f)Jjj の形に書けるということである.いずれにせよ パワースベクトル =L:ppに由来する部分 +L:qqに由来する部分 の形になる.右辺の 2 つの項は L 、ずれも周波数

f

の関数 で常に正の値を取る. Zpp に由来する部分 パワースベクトル

(

9

)

で定義される量をスベクトルに対する pp ブロックから のノイズ寄与率と名づける .L: の対角プロックの数が一 般の L 個であっても各ブロックからの寄与率を同じよう に定義することができる.

2

.

3

TIMSAC

TIMSAC パッケージに含まれている AR モデルにも とづく多次元時系列解析用のプログラムとして.

FPEC.

MULNOS.

DECONV 等がある. FPEC を使うことに

5

4

8

よって AR モデル (1)のパラメータの推定値が得られる. MULNOSは， FPEC の結果にもとづいてノイズ寄与率 (9)を推定するプログラムである.ただし， MULNOS は L: (4)が完全に対角型であること，つまり各雑音源が互い に独立であることを前提としている. DECONV は，雑音源から各測定値へのインパルス応 答を推定するプログラムであるが，本稿ではこれ以上触 れる余裕がない.

2

.

4

ARDOCK

ARDOCKは. MULNOS に手を加えて，一般の対角 型の分散共分散行列 L: (4)に対するノイズ寄与率(9)が推定 できるように改良したプログラムである. ARDOCK の機能を以下のように整理することができ る. [ARDOCK の繊能1 ノイズの分散共分散行現u: のプロヴク対角化:r;の構 造を変えたときの A

1

C(5)の値が計算される .L: の構造 が変わることによってパワースベクトルの推定値，ノイ ズ寄与率が変わる.

A 1

C の値が小さい場合のものが最 も信頼できるはずだが，問題意識のあり方によっては， AIC 最小のモデルを使えない場合もある.使用者の判 断で処理しなければならない. マスキング操作: 1. 周波数応答関数を計算する (6)式における AR モデ ルの係数行列 {A(m)} の任意の成分を強制的に 0 に置く ことが可能.たとえば，マスク 1110000 1110000 1110000 1111111 1111111 1111111 1111111 を (7x 7)係数行列 {A(m)} にかけることによって，

Ajk(m)=O

(j

=1.2

,

3

,

k=I.2

,

3

.

4

,

m=I

,

2.

…

.M)

とした場合を調べることができる.このような机上実験 によって，ある信号経路を切り離した時のシステムの挙 動の変化を容易に確認することができる.

2

.

L: の任意の対角プロックを強制的に O にできる. これによって，あるノイズ源、が除去された時の効果を確 認できることになる. パワースペクトル表示: 各変数のパワースペクトルのグラフを見ることができ る.細か L 、ことだが，対数スベクトルではなく， r 生の」 パワースベクトルである.ノイズ寄与はパワースベクト

ji;川川山

i

j::山川山門

l

iJN川い4"VI

出 jj山川附酬川川

3::仙川凡川川作オ

0::

4

1

V •

¥

J

¥

f

t

^^"'I"<品J 可If

¥

.rJ'Iv v ..

V

I.f\,.../

υi r lI¥

,

.rv

"J " V \.r ・""，-..

;::ド刊「γ

f二

Time 図 1 ロータリーキルンの運転記録 ルカ丸、くつかのCiE) 項の和であると L 、う見方に立つも のだから，これとを並べるパワースベクトルは「生J の ものであるべきなのである.グラフ表示に当っては，ス ケーリングを固定して，モデル操作にともなうスベグト ルの形の変化だけでなく，高さの変化がすぐわかるよう になっている.

3.

数値例

図 1 はセメントのロータリーキルンの運転記録であ る.長さ 511 ， 7 チャンネルである.第 1 変数から第 3 1989 年 10 月号 変数までが被制御変数，つまり被制御システムであるセ メントのロータリーキルンの状態を表わす変数，残りの 4 変数が制御入力である.図中 Wattage などと記され ている略称の意味に関しては表 1 にまとめておいた.な お， グラフに入れる目盛りの数値の桁数を小さくするた めに源データを 100 で割ってある.

[MODEL 1

]

AR モデルを当てはめたところ，次数 4 の AR モデルが AIC の意味で最良のモデルとして選 ばれた.ノイズ源の分散共分散行列を正規化した相関行 列は表 2 のようなものであった.非対角項が比較的小さ (37)

5

4

9

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

相関行列 表 2 数 変 表 1 1.00 -0.15 ー O.16 O.10 ー 0.06 0.15 -0.03 -0. 15 1.00 ー 0.01 0.08 0.01 0.03 -0.09 -0.16 ー 0.01 1. 00 0.04 0.09 0.02 -0.03

O

.

10 0.08 0.04 1.00 ー 0.11 0.09 -0.01 -0.06 0.01 0.09 ー O.11 1.00 ー 0.05 -0.06 0.15 0.03 0.02 0.09 ー 0.05 1.00 0.04 -0.03 -0.09 -0.03 -0.01 -0.06 0.04 1. 00 容 キルン駆動所要動力 クーラ下室圧力 窯尻ガス温度 キルン回転速度 燃料供給率 クーラグレート速度 窯尻ダンパ開度 内 略称

Wattage

Pressure

Temperature

R

Fuel

CR

Damper

変数番号 'ι 内 49Jd せ ZJZO ヲ g 雑音によっている.低周波領域における第 1 変数と第 3 変数からの寄与が大きいのが目だつ. 第 3 変数のスペクトルの構造も第 1 変数のものと似て いる.パワーの大きな部分が被制御変数の発生する雑音 に由来するものである.直流成分における第 6 変数から の寄与が顕著である. 制御変数のパワースペクトルに関しては，第 6 変数と 第 7 変数が比較的被制御変数の雑音に応答している以 外，ほとんど自分自身で閉鎖的な挙動を示しているのが 印象的である.

[MODEL 3J

ノイズ寄与率の定義を拡張してあるの で，必ずしも 7 変数全部が互いに独立でない場合にもノ イズ寄与率を推定することができる.制御変数のノイズ と，被制御変数のノイズのあいだには相関がないが，被 制御変数のノイズ相互の間に，あるいは制御変数のノイ ズ相互の聞には相関があると想定するモデルを作ってみ よう. MODELl の分散共分散行列にマスク くノイズ源の問に顕著な相関は見られない.このモデル を MODEL 1 と名づけよう.

MODEL

1 の AIC の値 は表 3 の第 1 行自に示されている 25765.92であった.こ れより次数の高いモデルも，低いモデルも AIC の値は これより大きくなる.このモデルから推定される各変数 のパワ}スベクトルのグラフを図 2 の 1 段目に示す. パワースベクトルを見ると，第 1 変数から，第 3 変数 までの被制御変数，および制御変数である第 4 から第 7 変数までのいずれをとっても， ほとんど周波数 0.1 以下 の低周波部分にパワーが集中していることがわかる.ス ペグトルの構造は単純である.

[MODEL 2

J

このスベクトルの構成を調べるために， MODELl の分散共分散行列にマスク 1000000 0100000 0010000 0001000 0000100 0000010 0000001 1110000 1110000 1110000 0001111 0001111 0001111 0001111 をかければよい.相関行列が 3

x

3 と 4

x

4 のプロッグ 対角型になる.これを MODEL3 と名づける.このモデ ルの AIC の値は表 3 の 3 行目に示す25776.73 となる. まだ，全変数の聞に相関有りとする最初のモデルより値 が大きく，不満は残るが，全変数が完全に独立であると する 2 番目のモデルに比べるとだいぶ低くなる.図 3 の 1 段目にこのモデルから推定される最初の 3 変数のパワ 対角ブロック型モデルの AIC 表 3 25765.92 25799. 14 25776. 73

AIC

MODEL1

MODEL2

MODEL3

分散共分散行列のプロック構造 ( 1-7) (1)(2)(3)(4)(5)(6)(7)

(

1

-

3

)

(

4

-

7) モデル をかけて MODEL 2 を構成して調べてみる.このモデ ルから推定される各変数のパワースベクトルは図 2 の 2 段目に示されている. 1 段目の結果とほとんど変わらな いなかで，第 2 変数のスベクトルの形の変化がマスクの 影響によるものである. MODEL2 の AIC の値は，表 3 の 2 行自の 25799. 14で，情報量規準の立場からは，全 変数が互いに無相関と見ることに無理があることを示し ているが，ここではあえてこのモデルにもとづくノイズ 寄与率を図 2 の 3 段目に示す. 第 1 変数のスペクトルの構造を見ると，直流成分のほ ぼ50%が制御変数の揺らぎに起因するものの，パワーの 大部分が自分自身の発生する雑音によっていることがわ かる.制御変数の中では，第 5 変数と第 7 変数からの寄 与が大きい.第 4 変数と第 6 変数からの寄与は小さい. 第 2 変数のスベクトルの構造は第 1 変数のものと似て いる.やはりパワーの大きな部分が自分自身の発生する

5

5

0

Wattagc Prcssurc Tcmpcratul'C R

凸

C J

円

kEl 6

凸

0

2

J同

F 4 3 -;::;• 同 Q 0

品2

4 _dQ

。

₌

₌

₂ ~ 4 ~ 2 ) ) 5 2

世て〉ωロCコ

P'H >qU

m自司

2 _{4D 25 1}

tL。

Fh >o' 4

=注ι。。ぃ

4 EECh CLE J

仏。ι。

E-Ze 。 0.10 0.20 0.0 0.10 0.20 0.10 0.20 0.0 0.10 0.20

Freq Freq Freq Freq

Wattage Pressure Temperature R ( 小3

2同口。:雪

3 ( 小1

_N

JKJ 6

日凸

C 2 4

₎

首Q

₀

₂

₂

口。:;:

4 2

“

〉ぴ【ー

?

、

ぺe〉-

_{m z}

qコ

_d J

、

j 2 ¥

判~自q3

4 J

盲

1 でqコ 2

E三

tHCh ECE B tLcbv

」

o, 4 p 。 。 0.10 0.20 0.10 0.20 0.0 0.10 0.20 0.0 0.10 0.20

Freq Freq Freq Freq

0.10 0.20 Freq ハ U nF 臼 八 υ l 1 3 [ l リ目 。 F 。 10 0.20 1プ】 C'l 0.10 0.20 Freq 図 2 各変数のパワースベクトルとノイズ寄与率 ースペクトルのグラフを示す.この場合，パワースベク トルに対する各変数のノイズからの寄与を分離して推定 することはできないが， r被制御変数群J からの寄与と 「制御入力群」からの寄与は分離できる.図 3 の 2 段目 にこの場合のノイズ寄与率のグラフを示す.図 2 とは， 各ノイズ源からの寄与が細かく分離されていないだけで なく， r被制御変数群 j と「制御変数群J へのパワーの 振り分けに関しても多少違っているはずだが，被制御変 数と制御変数の聞の絡み具合いを見るにはこれで十分で あろう.

[MODEL 4]

今度はこのモデルをいじってみること にしよう.まず，制御変数のノイズを 0 にしてみる.分 散共分散行列にマスク 1989 年 10 月号 1110000 1110000 1110000 0000000 0000000 0000000 0000000 をかけて得られるモデルを MODEL 4 とする.このモ デルから推定されるパワースベクトルを図 3 の上から 3 段目に示す. 1 段目のスベクトルに比べて背が低くなる ことがわかる.特に直流成分に近い低周波成分における パワーの低下がL 、ちじるしい.突は，この結果はわかっ ていた.図 3 を見ると各パワースベクトルにおける制御 プロックからの影響の割合がわかるから，この分がなく なった時のパワースベクトルの姿は，予想できるのであ (39)

5

5

1

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

Fuel CR Damper 10 100 円 同 C E

d凸

10

₎

。

_J

KQ

=

=

J

5

会 5

S 5 4ZC05 P 500

&

_向~

ω

v

。

-Lp 4

¥

Eh CEOZ 4 J 。 0.10 0.20 0.0 0.10 0.20 0.0 0.10 0.20

Freq Freq Freq

Fue! CR Dampef 10 100

c

:

;

-

(

。ュ ) Q

。冨

t

;

l

10 k0

J

Q 2 l

3

円.叫g〉司切qコ

J

、

u

一

k

忍bω

50

且hG。

>F J 4

¥

ES

ω

B

」

-OE4

¥

ECE。同

H Z J 。 l 、』 0.10 0.20 0.0 0.10 0.20 0.0 0.10 0.20 Freq Frcq Freq Fue1 CR Damper Line Samp!e 1.0 1.0

'

"

( 小3 Q

。

J

K2 J

。'"口

₌

_J

₌

J

i

o

s

.D

同ロ

0.5 {¥

三・制。内同ロ

0.5 。。 0.0 0.10 0.20 0.0 0.10 0.20 0.0 0.10 0.20 0.0 0.10 0.20

Freq Freq Freq Freq

図 2 (つつ.き) る.これは，制御系が出しているノイズを抑えることが できれば，システム全体の揺らぎを小さくできるであろ うということを示唆している.

[MODEL 5

]

では，関 1 の運転記録での制御は，有 害無益と考えるべきなのだろうか? この点をチェック するには，制御入力から被制御システムへの信号の経路 を切断してみれば良い. こんどは，分散共分散行列と AR 係数行列に，それぞ れ，マスタ 111 0000 1110000 1110000 1110000 1110000 1110000 0001111 1111111 0001111 1111111 00011111111111 0001111 1111111

5

5

2

(40) をかければよい.このモデルを MODEL5 と名づける. このモデルにもとづくパワースベクトルの推定値を図 3 の 4 段目に示す. 3 段目に比べてスベクトルの背が高く なった.直流成分のごく近傍以外の周波数領域では 段目のスベクトルより債が大きくなっている. これらの結果から，制御が確かに有効に働いてはL 、る が，その代償として非常に低い周波数領域における揺ら ぎを持ち込んでいるらしいということがわかる. まさ に，赤池らが図 1 のデータから読み取ったことである. これにつづく赤池らによるセメントのロータリーキルン の計算機制御の詳細に関しては参考文献を参照されたい

[MODEL 8

J

最後に，被制御システムから制御入力 への信号経路を切った場合のパワースペクトルとノイズ 寄与率のグラフを図 3 の 5 段目と 6 段目に示しておく. オベレーションズ・リサーチ

Temperature c。 J

w

g

4 2 主 2

3

も 一一一}包 .oom0.100J150.205

Freq

仏

PressllrC' M J

旨 3

0

.

2

3

1

O~ ，....:::: I ~ 0.0 0.05 0.10 0.150.20 }i Freq Wattage aaτ 。，“ 同 J 凶凸()冨)主 E8 包〉宣伝 6 0 0.0 0.05 0.10 0.150.20 Freq -;;:; Wa抗age 守 Pressure

f1:￨;

Temperature 1'emperature 0 0.0 0.05 0.100.15 0.20 Freq 建計 J

4

::;::

者 2

5

O~ -ー一一一一-=l !l 0.0 0.050.10 0.150.20 b Freq Jl" Pressure 叶唱 J 同 日 。

.

2 z;> 'u;

S

nJ 、ー'=

ﾖ

:

O

0.

'

0

5 0.10 0:15 0.20~ F四屯 巳4 3 Wattage Foa 告の 'u 唱、同凶 (52) 診百 5 司旬。=。内山

国

8

2

3

2

Freq 仏 Temperature 0 0.0 0.050.10 0.150.20 Freq

'

"

J 凶 。 。 2 z;>

2

4

5

O~ 、』← 0.0 0.05 0.10 0.15 0.20

i

i

!

o Freq :c... 4 Pressure 3 2 Wa抗age noa 告の FU 的 J 凶口。 EE28EE 。仏 シミュレーション このプログラムを使用することによって，ノイズ寄与率 および，各コンポーネントの聞の接続を切った場合のス ペクトルの変化を見ることができる. ある信号経路を切ったときの効果が小さければ，そこ は本来切れていることが予想される.これを確認するこ とは制約された AR モデルのあてはめを行なって，その モデルの AIC と無制約 AR モデルの AIC を比較する ことによって可能であろう. ftIJ御系を別のものに切りか えたときの効果を見るようにすることも容易だが，ここ ではこの点には深入りしないこととする. このように，モデルのパラメータを変えてみた場合の 結果をただちに信用するのは危険であるが，さらに深く 追求して L 、く場合の手がかりになり停ることは確かであ (41)

5

5

3

図 3 分散共分散行列と AR 係数行列に，それぞれ，マスク 1111111 1111111 1111111 0001111 0001111 0001111 0001111 線形ジステムにおけるコンポーネント聞の関連を調べ ることを目的とするプログラム ARDOCK を紹介した. 1989 年 10 月号 をかけた場合で‘ある.オベレーターがシステムの状態を 見ないで制御入力を操作したとすればこうなるであろう と L 、う非現実的な例である.

後

1110000 1110000 1110000 0001111 0001111 0001111 0001111

最

4

.

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

T

e

m

p

e

r

a

t

u

r

e

<0 J ~

4

2

者 2

5

"J -一一一一-a E も 00.050100.150.205

F

r

e

q

.

P

r

e

s

s

u

r

e

<0

S

3

0

き 2

2

3

O~ ，..-:::::" I

2

0

.

0

0

.

0

5

0

.

1

0

0

.

1

5

0

.

2

0

1

5

F

r

e

q

p..

W

a

t

t

a

g

e

6 4 2 3 凶凸()苫 )hdE 与む出。仏 ハリ ワ b A U PHυ l n u

oq

T i u -1

0

F

F 、u ) ( -n リ ハ U

•

ハ UAU

公 Wa枕 age 公 Pressure っ

y

iJ

T

c

m

p

e

r

a

t

u

r

e

(つつ'き)

[2J

八木原彬股( 1976) ，セメントプロセスの制御，数

理科学，

No. 153

,

p

p

.

5

3

-

5

9

.

[3 J Otomo

,

T.

,

Nakagawa

,

T. and Akaike

,

H.

(1972)

,

S

t

a

t

i

s

t

i

c

a

l

approach t

o

computer con

t

r

o

l

o

f

cement rotary kilns

,

Automatica J

.

IF

AC

,

8

,

p

p

.

3

5

-

4

8

.

[4 J Hagimura

,

S.

,

Saitoh

,

T. and Yagihara

,

Y.

(1988)

,

Application o

f

time s

e

r

i

e

s

a

n

a

l

y

s

i

s

and modern c

o

n

t

r

o

l

theory t

o

the cement

plant.

,

AISM

,

Vo

l

.

40

,

No.3

,

p

p

.

4

1

9

-

4

3

8

.

[5 J Nakamura

,

H. and Toyota

,

Y. (1988)

,

Staｭ

t

i

s

t

i

c

a

l

i

d

e

n

t

i

f

i

c

a

t

i

o

n

and optimal c

o

n

t

r

o

l

o

f

thermal power plants

,

AISM

,

Vo

l

.

40

,

No.1

,

p

p

.

I

-

2

8

.

[6J

尾崎統編( 1988) ，時系列論，放送大学教育振興 会.

ﾗ

ﾗ

ﾗ

図 3 る. Cプログラムの公開について 1 汎用性のあるプログラ ムで使い道もあると思われるので，なるべく早い機会に 何らかの形で公開することにした L 、が，時聞がとれるか どうかが問題である.

5

5

4

(

4

2

)

ここで紹介した方法については， 1986年度の統計数理 研究所共同利用研究集会 (61 ・共会・ 18) で，ある程度 発表済みである.その後の明治大学の大矢多喜夫教授と の「雑談」の中で刺激されて，プログラム化した部分も 多い.記して感謝の意を表したい. 辞] 文 [ 1 ] 赤池弘次，中川東一郎:ダイナミツクシステムの 統計的解析と制御，サイエンス社 (1972).

ﾗ

献

ﾗ

ラ老 [謝 参