時系列の三角関数分解 -- 最小2乗法

時系列を三角関数に分解するには、 三角関数の係数を最小2乗法によって求める方法もある。 これは、現象の周期はよくわかっているが、 観測の時間間隔が一定でなかったり、 観測が欠けていたり、まちがいがあると判断して除外する値があったりする 場合に適している。

ここでは、1年間の毎日の気温のデータをもとに、1年周期から1/6年周期 までの振幅と位相を決めたプログラム(Fortran)を紹介する。

• tsinfit.f がメインプログラムである。データを読みこみ、 次に説明するサブルーチンsinfitを呼び出し、 得られた結果を書き出す。 データを読むところは入力データの形式に合わせてなおす必要がある (ここでは下に述べるテストデータに合わせてある)。
• sinfit.f が最小2乗法を実行している主要部分である。 時系列の時間平均値ave、分散var、各周期成分の振幅amp、位相pha、 各周期成分の分散の全分散に対する百分率prcvarを計算する。 なお、計算不可能な場合は、これらの変数にはqndとして指定された値がはいる。 線形最小2乗法のプログラムとしてはなるべく一般的に作ったつもりだが、 次の点がこの問題専用なので、他の問題に応用するにはなおす必要があるだろう。
  • 配列tの寸法nmax (n以上であること)
  • 位相を、ここでは、温度が最高となるような1月1日からの通しの日づけで 表現しているが、一般には角度(度あるいはラジアン)で表現する。 dpha = dpha * tcycle/k / (2*pi)を 取り除けば、ラジアンのままになる。
  • 次に述べる連立方程式が解けるかどうかのしきい値epsとして、10^-6という値を設定しているが、これがいつも適切かどうかは確かめていない。
• gauss.fは、 連立1次方程式を解くサブルーチンである。 これはGaussの消去法(LU分解、後退代入)と呼ばれる方法で、 村田(1986)の教科書の例文ほぼそのままである。 (ただし、変数を全部宣言し、変数名を小文字で書いた。 また、DO文の表現を、Fortran 77標準の表現から、現在多くのコンパイラで利用できる END DOを使った表現に変えた。) 解くべき方程式の係数行列が対称であることを考えると、 コレスキー(Cholesky)法のほうが適切だが、 計算量が特別多くはないので変えなかった。
• soya1997.txtは テスト用に用意したデータで、 気象庁アメダスの1997年の毎日の気象観測値 (気象業務支援センターから発行されているCD-ROMの内容を増田が形式変換したもの) のうちから、宗谷岬の観測値をぬき出したものである。 各列左端から51--56バイトめに、日平均気温が0.1度C単位で書かれている。 例文ではこの数値を読む。なお、32767は観測データが欠けていることを示す。
• tsinfit_sh は、このプログラムをコンパイルし、テストデータを読みこんで実行する手順の シェルスクリプトである。

最小2乗法の入門的解説を書いた本はたくさんあるが、 とりあえず、粟屋 (1983)を紹介しておく。 この本の副題「アナログとディジタル」は、 コンピュータ上のデータ形式ではなく、対象が離散量か連続量かという意味である。 この本のねらいは、放射線実験などのポアソン分布に従うと想定される離散量 を扱うために、最小2乗法に代わる方法を提案することにあるのだが、 最小2乗法の理解のためには、 この本の「アナログデータ」に関する部分だけを読めばよい。 三角関数をあてはめる例はないが、多項式の場合と同様に考えればよい。

最小2乗法を本格的に使うプログラムを書くのには、 中川・小柳(1982)が参考になる。 ただし、線形代数、統計学、数値計算アルゴリズムの基礎知識が必要 (わたしも全部を理解しているわけではない)。 また、この本は著者たちが中心となって作ったSALSというプログラムを 紹介しているが、これはかつての大型計算機の環境に合わせて作られたため、 今となってはあまり使われなくなってしまったようである。

参考文献

• 粟屋 隆, 1983: データ解析--アナログとディジタル。 学会出版センター, 207 pp. (その後第2版が出ている)
• 中川 徹、小柳 義夫, 1982: 最小二乗法による実験データ解析--プログラムSALS。東京大学出版会, 206 pp.
• 村田 健郎, 1986: 線形代数と線形計算法序説。サイエンス社, 225 pp.

[目次へ]