ブロック図+実例からARモデルを理解する

#はじめに

制御工学における大きなテーマの一つに、システム同定があります。

システム同定（システムどうてい、System Identification）とは、計測データから動的モデルを構築するための数学的ツールやアルゴリズムを指す用語

「動的モデル」の代表例としてARXモデル（Auto-Regressive with eXogenous model、自己回帰＋外乱外性入力モデル）が挙げられます。ARXモデルは専門書の冒頭に書いてあるようなド定番と言えます。

一方でARモデル(Auto-Regressive model、自己回帰モデル)というのもあります。ARXモデルからXを取ったものなので、これもシステム同定のためのモデルかと思い調べてみたらどうも違うらしい。統計学の時系列解析のためのモデルとのこと。

パラメータ同定に使うARモデル、システム同定とか信号処理とか統計学とか色々なところに出てきて結局出自はどこなのよ？　と思ったら足立先生のシステム同定本に「統計学の時系列解析で用いられる」「システム同定問題ではない」と書いてあって流石神ですねと思った。https://t.co/0xwW6U13rf

— モータ制御マン (@motorcontrolman) November 1, 2019

この流れでARモデル気になるな~と思い、真面目に勉強してみると「個人的に気持ち悪い部分」があったため理解に少し苦しみました。

本投稿では気持ち悪さの原因について述べ、次にブロック図をもとに気持ち悪さを解消、ARモデルを理解するまでを示します。最後にARモデルの実例を挙げ、なぜARモデルが使われているのか、どう当てはめているのかを考察することで理解を深めるものとします。

あくまで個人的解釈として書いているので、間違いがあればどしどしご指摘下さい。

#1.ARXモデル
ARモデルを説明する前に、ARXモデルについて説明します。式は下記。

A(q)y(k)=B(q)u(k)+w(k)　\\
\:\\
A(q) = 1+a_1q^{-1}+a_2q^{-2}+...\\
B(q) = 1+b_1q^{-1}+b_2q^{-2}+...

ここで、yは出力、uは入力、wは誤差やノイズを意味し、kは離散時間系におけるサンプル点を意味します。
qはシフトオペレータと呼ばれるもので、z変換における演算子zと意味は同じです。
ARXモデルのブロック図は下記です。

qをzに書き換えた上で下記のように表現すれば、制御屋にとってはなじみ深い図になるかと思われます。FBループとコントローラは本来のARXモデルにはないですが、書いた方がしっくりくるかと。

（追記注：あくまでのイメージであって、実際のモデル同定では観測雑音の影響を回避するためループを切って行う必要があります。Thanks to @wcrvtさん）

すなわちARXモデルは、コントローラによって決定される操作量とプラントからの出力を利用したシステム同定であり、ノイズは正確な同定を阻害する因子と言えます。非常に分かり易いですね。

2.ARモデル

一方のARモデルですが、ARXモデルのB(q)u(q)項をなくすことで得られます。

A(q)y(k)=w(k)

ARXモデルより簡単ですね、別に気持ち悪い感じもしません。
次に、ARXモデルの図を流用したARモデルの図を示します。

この図、どうでしょうか。冒頭で述べた「気持ち悪さ」はこの図に対して感じました。

この図からするとARモデルはノイズとプラント出力を利用したシステム同定に見えます。ARモデルに対する私の第一印象は、「操作量を使わず、ノイズ＝素性の分からない信号からシステム同定するなんて意味が分からない…」と言った感じでした。

3.ブロック図を変形、気持ち悪さを解消する

A(q)の定義より、ARモデルの式は下記変形できます。

(1+a_1q^{-1}+a_2q^{-2}+...)y(k) = w(k)\\
\begin{align}
y(k) &= (-a_1q^{-1}-a_2q^{-2}-...)y(k)+w(k)\\
     &= (-a_1-a_2q^{-1}-...)y(k-1)+w(k)\\
     &= A_1(q)y(k-1)+w(k)
\end{align}

A_1(q) = -a_1-a_2q^{-1}-...

図で表すと下記です。便宜上、qをzに置き換えています。

離散時間状態方程式のブロック図表現のような形が得られました。注意しなければならない点として、離散時間状態方程式のブロック図におけるAは行列ですが上図における$ A_1 $は伝達関数です。また離散時間状態方程式は他入力多出力を扱うことが出来ますが、ARモデルは単出力のみしか扱えません。

上記はARモデルのみの図であって、システム同定に関するブロックは含まれていません。システム同定は出力yの推定値と真値との誤差を利用するため、下記のように表せます。ここで$ \hat{A_1}(z) $は推定のための伝達関数です。

上図よりARモデルによるシステム同定は、プラント出力の現在値および過去値を利用した同定法であって、推定誤差が$ \hat{A_1}(z) $によって変化する関係を用いることで同定を行うものと言えます。言うまでもなくノイズwはシステム同定に用いられません。
(ちなみにARモデルでは$ \hat{A_1}(z) $の決定法として最小二乗法や最尤法などが挙げられます、本稿では省略。最小二乗法の考え方は線形回帰の場合とほぼ同じです)

上記の説明を見て「そんな変形しなくても分かったことでは？」という方もいらっしゃるとは思いますが、少なくとも私(制御屋)にとってはこのブロック図を書くことでやっとARモデルを理解できました。恐らくですが、「モデル」という言葉や元のブロック図に引っ張られすぎたのではと思います。

なお、統計・確率の分野におけるARモデルによるシステム同定のための式は下記です。この式が苦労なく理解できる方であれば、上記のブロック図は不要なものと言えるでしょう。

y(k)=\sum_{i=1}^{p} a_i y(k-i) +w(k)\\
\hat{y}(k)=\sum_{i=1}^{p} \hat{a}_i y(k-i)\\
\epsilon(k)=y(k)-\hat{y}(k)

さらにおまけとして、上記のブロック図は下記のようにも表せます。言わずもながら、IIRの形ですね。信号処理に親しい人であれば、この図のほうが自然に見えるかも知れません。

ARモデルに限った話ではありませんが、式であっても図であっても色々な表現方法があるので自分にとって一番分かり易いものを探したり、作ったりするのは深い理解の助けになります。この記事はその１例とも言えます。

ところで、ARモデルにノイズが無かったらどうなる？

重要なことなので、実例紹介に入る前に記載します。

前述の通りARモデルによるシステム同定は、プラント出力の現在値および過去値を利用したものであってノイズは同定に用いられません。

そうなるとノイズの存在は同定を邪魔するものに思われます。いっそのことノイズがないほうが正しい同定が行われるのでは？と思えます。

では、ノイズが無かったらどうなるでしょうか。
１例として$ a_1=-0.9 $ のARモデルの出力yおよびノイズwを下記に示します。

MATLABコードは下記。ツールボックスは一切不要です。

w = randn(1,100);
y = zeros(1,100);
a1 = -0.9;

for i=2:100
   y(i)=a1*y(i-1)+w(i);
end
plot(w);hold on;plot(y);
legend('w','y');

ノイズwに0をかけて、ノイズを無くしてみましょう。するとどうなるか。

当然のことながら、yはずっと0になります。これではシステム同定出来ません。当たり前といえばそれまでですが、これはARモデルにおいて非常に大事です。

すなわち、ARモデルはノイズ入力があるからこそシステム同定が可能なのです。ノイズというより、観測者が知りようのないランダムな入力と捉えるとよいかも知れません。（ここを理解しておくとこの後の事例が非常に分かり易くなります）

なおARXモデルは入力にノイズw以外の信号uを持つため、ノイズが無くてもシステム同定を行うことが可能です。ARXモデルにとってノイズは毒ですが、ARモデルにとってノイズは毒でも薬でもあります（ドヤ顔）

（追記注：PE(Persistently Exciting)性の観点からすると、ARXモデルにとってもノイズが薬たりうる場合あり。 Thanks to @wcrvtさん、コメント欄参照）

4.ARモデルの実例

①ファイナンス

ファイナンスのための数学基礎 第1回 オリエンテーション、ベクトル

上記は慶応義塾大学の長倉先生の講演資料ですが、ARモデル+経済学、ARモデル+金融工学で調べると数えきれないほどの文献が出てきます。恥ずかしい話ですが全然知りませんでした。

前述のようにARモデルは、観測者が知りようのないランダムな入力を取るシステムを同定するものです。ファイナンスにおける出力はお金の動きなので観測可能ですが、入力は色々な要素がありすぎて捕捉のしようがありません。なので、入力＝ランダムとみなしてARモデルを適用することでモデルの同定を行い、将来予想に用いているのだと思います。（間違ってたらごめんなさい）

制御屋の私にとって、システム同定のモデルだと勝手に思っていたARモデルがこれらの分野で使われているのは非常に興味深いです。また同時に、システム同定・時系列解析について学べばこれら分野で使われる数学も理解できるようになれると思うと楽しくなってきます。

②音声

音声の特徴抽出(ＤＦＴ, ＬＰＣ, ケプストラム分析)

上記は東京大学 高道先生のNAISTでの講義資料。ARモデルについてはp48以降を参照、LPC(Linear Predictive Coding)分析法として挙げられています。

音声分析にARモデルをあてはめる上での個人的に面白いところは、「音声は発生源においては様々な周波数成分を持つランダムな信号である」と捉えられることです。なぜかというと、ARモデルのA(q)に該当する部分は下記図における声帯以降と考えられるからです。(上記リンクのp24より引用)

人間が任意の声を出そうとしたときは、声帯以降の動きを変化させることでARXモデルのA(q)に変化が引き起こされ、結果として望み通りの声が出ます。

ではARモデルにおける「ノイズ」はどこから来るかというと声帯より前の部分と考えられます。この「ノイズ」が様々な周波数成分を持っているからこそ人間は色々な音声を出すことが出来る、ARモデルからはそのように考えられます。もしも様々な周波数成分を持っていなかったなら、人間が放つ音声は全然違うものになっていたと考えると興味深いです。

ちなみに著者は、音声処理にARモデルが使われているという知見はMathworksセミナーにおける下記のデモから間接的に知りました。ツールボックスにSignal Processing、Statistics & Machine Learning が必要な点がちょっとあれですが、非常に面白いデモなので是非お試し下さい。

マイクから読み取った音声の日本語母音認識

上記デモはプログラムの内にてユール・ウォーカー法によるパワースペクトル密度の推定を利用しています。ユール・ウォーカー法はARモデルにおける$ \hat{A_1}(z) $の決定法の１つです。

③ストレスの数値化

・ARモデルのパラメータ推定法 自己回帰モデル（AR法）を用いたパワースペクトル密度算出　その１
・自己回帰モデルより推定した応答関数による 心拍変動

心拍からの特徴量分析にARモデルを用います。特徴量はフーリエ変換によっても求められますが、ARモデルを用いるほうが利用しやすいためARモデルが広く用いられるとのこと。

ストレスの数値化にARモデルをあてはめる上での個人的に面白いところは、心拍はストレスによって直接変化するのではなくARモデルによって間接的に変化すると捉えられること、およびストレスによって変化するのはARモデルのA(q)であると捉えられることです。
分かりにくいので図にするとこうです。

音声と同様、心拍も元となっているのはランダムな信号なのです。(実態は分からないが、ARモデルを当てはめているのでそう取れる)
心拍がランダムってなんだかちょっと怖いですね。またA(q)に当たる部分が無かったら人間の心拍はランダムのまま出力されていたことでしょう。あって良かった、A(q)。

#おわりに
当初、ARXモデルを主体とするシステム同定について勉強するつもりだったのですがものすごい勢いでARモデルに脱線しました。
そこで感じた気持ち悪さについて考察することで、ARモデルが何なのか腹落ちすると同時に制御と時系列解析の距離感がなんとなく分かりました。距離感が分かったおかげで、時系列解析の文献を読む上でのポイントが分かったのが個人的には一番の成果でした。

本投稿を元に、読者が著者の感じたことを追体験し「ARモデルだいたいわかった」「時系列解析ってそうなのね」となって頂ければイイナ！と思います。