信号解析第２回講義録日時：２０１５年４月２７日（月）講義内容：時系列信号の確率・統計処理担当者：情報知能工学科小島史男 1 はじめに前回の講義では、信号に含まれる不確定性について説明しました。時系列信号の取り扱いの基本は確率・統計量です。今週の講義では、時系列信号の特徴を表す平均値関数、自己共分散関数、自己相関関数について説明します。酔歩というモデルについて先週説明をしましたが、まず確率・統計の基本をおさらいし、弱定常過程の時系列信号について考察していきます。 2 連続確率変数われわれの日常生活におけるデータを確率変数としてとらえたとき、どのような取り扱いが必要でしょうか。株価変動を時系列信号の例にとると、この変数の適用範囲は円という貨幣単位です。また１日の気温の変化の記録はどうでしょうか。温度の物理量は実数ということになります。また交流の電気信号なら複素数となるでしょう。ここではもっとも一般的な確率変数が Y ∈ (−∞, +∞) で値域が与えられる場合について考えましょう。この場合確率は Prob[Y ≤ y] = G(y) で与えます。ここで G(y) は以下の性質をもつ単調増加関数です。 0 ≤ G(y) ≤ 1, G(+∞) = 1, If y1 ≤ y2 , Then G(y1 ) ≤ G(y2 ) またこの関数が y について連続で微分可能とすれば次の関数を定義することができます。 g(y) := dG(y) dy これらの関数 G(y), g(y) をそれぞれ分布関数、確率密度関数と呼びます。これらを用いると酔歩と同様に、連続確率変数の平均値および分散はそれぞれ次のように記述することができます。 ∫ µ = E(Y ) = ∫ +∞ −∞ ydG(y) = Var = E{(Y − µ)2 } = ∫ +∞ −∞ +∞ −∞ yg(y)dy (y − µ)2 g(y)dy 確率密度関数としてよく知られているのは正規型ですね。これは { 1 (y − µ)2 g(y) = √ exp − 2σ 2 2πσ } で与えられ、平均 µ および分散 σ 2 という２つの母数（パラメータ）によって分布が決定されることはこれまで何度も学習してきたはずです。 3 時系列の統計量時系列 y[n] が実数値の確率変数で与えられるとき、前節の結果でその統計量を記述することができます。平均値や分散はそれぞれ µ[n] = E(y[n]), Var[n] = E{(y[n] − µ[n])2 } (n = 0, 1, 2, · · ·) と記述できます。これらは時間の関数 n に依存して決まりますので、それぞれ平均値関数、分散関数と呼ばれます。ところで、時系列では時間の経過のあいだの関係性を記述できることが本質的です。時間 n ごとの統計量だけでは、時間変化の統計的記述はできていません。これに関連する統計量として以下に示す自己共分散関数があります。 Cov(n, n − k) = E{(y[n] − µ[n])(y[n − k] − µ[n − k])} 時系列の統計量としては今後、平均値関数、分散関数に加え、自己共分散関数という特徴量を考えていくことになります。 4 弱定常過程 −1 いまデータ長 N の時系列信号 {y[n]}N n=0 が得られたとします。このときこの時系列の統計量としては、 {y[0], y[1], · · · , y[N − 1]} が同時に起こる同時分布関数、同時確率密度分布を考えないといけないことになります。しかしこれは問題を大変複雑化してしまうので、すこし制限を設けましょう。いま、時 −1 系列 {y[n]}N 、そして自己共分散 n=0 の平均値関数、分散関数が時間 n に依存しない（すなわち定数値）関数は時間差のみに依存するとします。すなわち、任意の整数 l に対して µ = E(y[n]) = E(y[n − l]) C0 = Var[n] = Var[n − l] Ck = Cov(n.n − k) = Cov(n − l, n − l − k) が満足される時系列信号を弱定常確率過程と呼びます。さらに自己共分散関数を規格化した自己相関関数 Rk = Ck C0 は大変便利な指標です。この値はかならず −1 と +1 の間で値をとります。正数の場合“ 正の相関 ”、負数の場合“ 負の相関 ”、また絶対値が 0 に近いほど関係が薄く、1 に近いほど関係が強いという評価ができます。 5 時系列データの統計量の推定では本日の講義でもっとも重要な部分を学習します。これまでの講義で時系列信号の統計量について、それがもし弱定常過程ならば平均値（定数）、自己共分散関数（時間差の関数）、自己相関関数 (値域 [−1, 1]) で記述できることがわかりました。でも少し考えてみると、われわれが日常取得できる時系列がはたして弱定常過程なのでしょうか。これはよく考えていかないといけないことですがこの問題はこれから少しずつ考えていくことにします。当面われわれが取り扱う時系列信号は実数値弱定常過程としましょう。それでもまだ問題があります。平均値、自己共分散関数、自己相関関数を計算するには、結局確率密度関数 g(y) を知っている必要があります。しかしわれわれのまわりの信号の分布はわかりませんので、何らかの方法で結局推定せざるを得ません。そのひとつの方法は、分布の形を仮定することです。確率モデルとはそのようなもので、たとえばガウス型時系列分布としてのブラウン運動はその典型例です。もうひとつの方法は時系列データから直接推定する方法です。実数値をとる定常時系列信号 −1 {y[n]}N n=0 の推定統計量として以下の３つの特徴量を紹介します。 (1) 平均の推定値（標本平均） µ ˆ= −1 1 N∑ y[n] N n=0 (2) 自己共分散関数の推定値（標本自己共分散関数） −1 1 N∑ Cˆk = (y[n] − µ ˆ)(y[n − k] − µ ˆ) N n=0 (3) 自己相関関数の推定値（標本自己相関関数） ˆ ˆ k = Ck R Cˆ0 ˆ k は時系列信号から直接計算できるので、それぞれの統計量に標本というこれらの推定統計量 µ ˆ, Cˆk , R 接頭語をつけて、標本平均、標本自己共分散関数、標本自己相関関数と呼びます。これらの推定量はサンプリングデータから直接計算できますが、これがなぜ意味をもつのでしょうか。教科書にはあまり詳しく書いていませんが、弱定常確率過程にエルゴード性を仮定すれば、この標本数 N を無限にとっていけば、前節で紹介した統計量に確率 1 で収束することが保証されています。したがって弱定常過程であれば、上記３つの統計量が十分よい近似を与えることになります。どちらにしても時系列データから直接計算できるのはこの３つの統計量ということに注意してください。 R をつかっておさらい 6 時系列信号の 4 つの標本統計量は R においてはそれぞれの組み込み関数 mean(y),var(y),acf(y, type = ”covariance”) および acf(y, type = ”correlation”) を使って求めることができます。教科書 p41 の式 (3.18) の信号生成モデルにもとづきこれまでの復習をしましょう。時系列信号が ( y[n] = cos 2πn 10 ) ( + cos 2πn 4 ) + w[n] で与えられるとします。データ長は N = 400、また w[n] は平均値 0、分散が正定値の正規型白色雑音とします。これは２つの周波数成分を伝送したときに途中で雑音が混入した状況を表しています。R では次の関数で記述することができます。時系列信号のモデル関数 model <- function(n,var) { y <- numeric(n+1) w <- rnorm(n, sd=sqrt(var)) pi <- acos(-1) y[0] = 2 for (i in 1:n) { y[i] = cos(2*pi*i/10) + cos(2*pi*i/4) + w[i] } return(y) } この関数は引数として、データ長 N = n、また白色雑音の分散 var を指定します。関数の中身ですが、白色雑音 w[n] の生成は、組み込み関数 norm(n, mean, sd) を使っています。これは正規乱数の生成関数です。先週酔歩では一様乱数 runif() を使いました。このように R では乱数の生成関数は接頭語として ”r” を、そのあとに確率分布の略号が続きます。ここではそれぞれの時系列信号のシミュレーションおよびその標本自己相関係数（自己相関係数の推定値）を表示しましょう。R でのコマンドラインは次のようになります。コマンドライン >source("C:/model.R") >y <- model(400,0.01) # for Figure 1 # >y <- model(400, 1) for Figure 2 >par(mfrow=c(2,1)) >plot(y,type="l") >acf(y,type="correlation") >dev.copy2eps(file="model-1.eps") # >dev.copy2eps(file="model-2.eps") 0 −2 −1 y 1 2 命令文のなかで acf が自己相関係数です。補足ですがコマンドライン中の par(mfrow = c(2, 1)) はグラフを同時に２個縦に並べるという意味です。また関数 acf はデフォルトでは相関係数（相関関数の場合は type = ”covariance” と指定）でグラフ表示が自動的に実行されます。これを実行した結果の例をみましょう。図 1 は白色雑音の分散が 0.01、図 2 は白色雑音の分散が 1 の場合の結果です。 0 100 200 300 400 Index −0.5 ACF 0.5 1.0 Series y 0 5 10 15 20 25 Lag 図 1: 時系列信号の標本過程とその統計量 (1) みなさんには次の宿題をだしますので、各自確かめてみてください。レポートはＡ４用紙に記入してください。提出先は自然科学系先端融合研究環（総合研究棟）３号館３０１号室です。部屋の前に回収ボックスを置いておきますのでそこに入れてください。また質問がある場合は、あらかじめ、[email protected] まで問い合わせください。なお R の実行ファイルは 2 0 −4 −2 y 0 100 200 300 400 Index 0.5 0.0 −0.5 ACF 1.0 Series y 0 5 10 15 20 25 Lag 図 2: 時系列信号の標本過程とその統計量 (2) http://cran.r-project.org/ からダウンロード可能です。宿題： (1) 英国のガス使用量の時系列データ UKgas を用いてその標本統計量を調べてみましょう。R において y < −UKgas としてあとは上記と同じ手続きを実行すればできます。 (2) 標準正規型白色雑音 µ = 0, σ 2 = 1 の疑似乱数を時系列（y < −rnorm(n, sd = 1), n は適当なデータ長）としたときその標本統計量を示してください。これは (1) と比べてどのような特徴があるでしょうか。