信号解析　第６回講義録

日時：２００６年５月２２日

講義内容：最小自乗法とデータ曲線

担当者：情報知能工学科　小島史男

1. はじめに

　先週はカルバックの情報量について説明し、データから真の統計モデルを推定する評価法について学習しました。今週はまず統計モデルの構造を与え、それに含まれるパラメータをデータから推定する手法について説明します。受講生のみなさんは、最小自乗法をよく知っているとおもいます。今週はまず最小自乗法が統計モデルとしてガウス型を仮定したときの最尤推定と一致することを示し、次にデータから最小自乗法により回帰曲線を求める計算法について説明します。赤池情報量基準(AIC)は回帰曲線を求めるメタ最適化手法として有名です。開学記念日で１回抜けましたので、今週で信号解析前半が終了します。説明する事項がたくさんありますが、基本的事項はしっかりおさえていきましょう。

2. 最大尤度法とパラメータ推定

　前回の講義で、尤度（対数尤度）は時系列信号の背景にある統計量を評価するもっとも重要な指標であることを説明しました。尤度関数はよく知られているこれまでのいろいろな計算法に明確な解釈を与えることができます。たとえば最小自乗法は誰でも知っている統計処理法ですが、その平均の意味を尤度関数から考えてみましょう。対数尤度を計算するには、統計モデルの構造を仮定して、その関数の係数をパラメータとしてこの値を変化させることにより、いろいろな尤度の値を調べる方法が一般的に用いられます。すなわち対数尤度をの関数

として表記します。この統計モデルをパラメトリックモデルと呼び、対数尤度を最大にするパラメータを求める計算法を最尤法(maximum likelihood method)と呼びます。では以下で統計モデルが平均値分散１のガウス型確率密度関数

で与えられるときに、平均値を最尤推定により計算してみましょう。この場合の対数尤度は前回の講義でおこなった計算を施すと

となります。この対数尤度を最大とするは、上式において、に関する極値をとり、

によって計算できます。時系列信号が独立に得られる場合は、データ長を無限大にとれば、その標本平均の極限の分布が平均値分散１のガウス型の統計モデルに一致することになります。このパラメータ最適化（最大尤度法）は自乗誤差

の最小化にほかなりません。これはよく知られている最小自乗法(least square method)ですが、その最適値である標本平均は以上の意味で合理的に説明できます。最尤法に関する計算アルゴリズムの詳細は講義の目的ではありませんので、ここではこれ以上深入りしないことにします。

３. 回帰モデルによるデータのあてはめ

　説明に入る前に、まず例題をみてみましょう。次の図はいま話題の阪神電鉄の株価１年間の毎日の最終終値を時系列的にまとめたものです。青色の不規則な変動が実際のデータ、ピンク色の曲線が最小自乗法により求めた回帰モデルによる曲線です。

このようなことはどうしてできるのかが今回の講義の目的です。ではまず回帰モデルについて説明していきます。いまデータの変動をとあらわし、各時刻におけるデータが個の変数の線形結合と不確定要素との和で与えられるとします。すなわち

のことを回帰モデル(regression model)と呼びます。回帰モデルではを目的変数、をこのデータが生成する説明変数と呼びます。説明変数の個数は次数(order)と呼び、説明変数の前についている係数を回帰係数(regression coefficients)と呼びます。第２項のは説明変数では記述できない不規則変動をモデル化したもので、残差(residual)と呼びます。ここで残差の統計モデルは平均分散の独立なガウス分布で与えたとします。以上の設定のもとで、適当な次元と説明変数を具体的に与え、パラメータベクトル

を適当に与えれば回帰モデルが計算できます。残差が互いに独立な分布であることを考慮すれば、尤度および対数尤度はそれぞれ

と表記できます。また統計モデルを平均分散のガウス分布としたので、尤度関数、対数尤度関数の確率密度は

となります。結局対数尤度関数

が導かれます。この尤度関数を最大にするパラメータベクトルを求めることにより統計モデルのもとで最良な近似値を求めることができます。たとえば分散に関して、上記の対数尤度を偏微分してみましょう。

この極値を求めると、

となりますから、各時刻の時系列データをし、任意の回帰係数を与えれば分散の推定値が求まります。一方を対数尤度関数に代入すると、対数尤度は

で与えられますので、これをについて最大にするには、分散の最尤推定量を最小にすればよいことになります。すなわち前節で述べたように最小自乗法そのものとなります。

４．ハウスホルダー変換による計算アルゴリズム

前節の結果から、回帰係数の最尤推定量は残差の最小自乗法解であることを示しました。ここでは、その解法アルゴリズムについて簡単に説明しましょう。まず回帰モデル全体を行列でまとめていきましょう。次元ベクトル、行列を

と表記すると、回帰モデルは簡単に

と記述できます。つまり最尤推定量はベクトルについて以下の次元ユークリッドノルムの最小解を求めることと同等です。

これを直接求めてもよいのですが、計算を効率よく実行するために直交行列を用いて以下のように変換します。

このような変換をしてもノルムの値は変化しませんから、変換後の最小化問題を解いてもよいことになります。この最小解はに関する線形代数方程式

を解けば良いことになります。を解きやすい形に変形できればアルゴリズムは簡単になります。そこで考えられたのが、ハウスホルダー変換(Householder transformation)です。

説明変数の行列の右側に目的変数ベクトルを加えた拡大行列を

のように定義します。この行列に適当なハウスホルダー変換を施して上三角行列

に行変換します。これを利用すると

と変形できて、この式から最小解は下記の線形代数方程式を解けば良いことになります。

後退代入をもちいることにより、回帰係数は後ろから順次求めていくことができます。以上の計算アルゴリズムをまとめてみましょう。

また残差の分散の推定量は

によって計算できることになります。

５. 赤池情報量基準(AIC)

前節までの議論で最初の例題を解くほとんどの準備が整いました。最後にひとつだけ問題があります。回帰モデルの次元はどう決めたらいいのでしょうか。データが多ければ大きくとればよいのでしょうか。もしこの考えでいくと、前節の行列のサイズがどんどん大きくなり、そのうちに実現不可能な大きさになってしまいます。できれば小さいほうが望ましいですね。このパラメータの次元決定に関しては、有名な赤池の情報量基準（Akaike Information Criteria, 以降AIC基準と呼ぶ）があります。この方法は実際には計算できない平均対数尤度と時系列データにより計算した対数尤度とのバイアス値（偏り）を最適化するために求められた指標です。具体的な計算は省略しますが、AIC基準とは

をモデル選択の基準として用いる方法です。ここで右辺最初の項は最大対数尤度ですが、２項目のkはパラメータの次元です。AICを最小にするモデルを選択することにより、できるだけ少ないパラメータ数で不規則信号の背景にある統計量を合理的に設計できる手法であり、簡便な手法であることから広く用いられています。ではここでの問題にAICを適用してみましょう。前節の計算アルゴリズムにより最大対数尤度は

であること、およびパラメータの次元がであることに注意すれば、この回帰モデルの赤池情報量基準は

となります。回帰モデルの次元は小さいほどよいのですから、実際は適当な最大次元（通常は20前後に設定）をさだめ、として１から順番に次元をあげていってが最小値となったところで計算を止め、そのときの次元を最適なとして、そのときのみ線形代数方程式を解き回帰モデルをつくります。

６．例題

　では最後に、今回の講義で最初に示した回帰曲線のあてはめを実行しましょう。ここでは参考図書の５章で実際に解いた例を参考にします。回帰モデルを

で与えます。ただしはまたはとします。この場合説明変数およびパラメータを

のように与えられることにします。また説明変数に含まれるですが、およそ1年の周期変動を考慮してと設定します。また考える最高次数とするとということになります。このようにして、次数を０から２１まであげていき、それぞれのAICを計算します。そのうちの最低の値をもつ次数のときの回帰モデルを最適として回帰曲線を当てはめます。表１はAICで求めた最適な次数です。また表２にここで求めた回帰モデルの詳細について掲載しておきます。最適な次数が考え得る上限値に達してしまいました。この場合は最高次数を以上あげていくと別の解がでてくるかもしれません。興味のある受講生はトライしてみてください。プログラムはダウンロードして使ってください。なおもし自分でプログラムをつくるのでしたら、ハウスホルダー変換の詳細については

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

技術評論社 ; ISBN: 4874085601 ; (1993/06)

を参考にしてください。

表１：最適な次数とAICの値

最適次数	AIC	推定分散
21	2492	1179

表２　回帰モデルの残差の分散およびAIC(m)の値

次数	残差分散	AIC	最小AIC値との差
0	685303	4022	1530
1	66540	3448	956
2	55347	3405	913
3	7620	2917	425
4	7199	2905	413
5	5553	2843	351
6	5312	2834	342
7	4790	2810	318
8	4262	2783	291
9	3887	2763	271
10	2960	2697	205
11	2959	2699	207
12	2847	2692	200
13	2069	2615	123
14	1997	2608	116
15	1891	2597	105
16	1886	2598	106
17	1566	2554	62
18	1326	2515	23
19	1325	2517	25
20	1324	2518	27
21	1179	2492	0

コメント：連続2回休講となりますので、今回はすこし無理して前に進めてしまいました。講義で十分説明できなかったことを講義録にアップしましたので、内容をよく確認してください。これで信号解析半分が終わりましたので、次回２９日は中間テストを実施します。時間は60分を考えています。ノート、ホームページの印刷は可としますので、じっくり考えて解答してください。また中間テスト問題の最後のところに次回講義（6月12日）までの宿題をだしておきますので、レポート作成してください。前半の講義で不規則時系列信号の基本的説明をおわりました。後半では時系列信号のモデルであるARMA過程を紹介し、そのあと状態空間モデルの推定問題としていまひろく使われているカルマンフィルタについて考えていきます。