Re:ゼロから始めるML生活

どちらかといえばエミリア派です

【R】時系列分析の覚書(導入)

こんな本を買ってみたので、エンジニアというよりデータサイエンティストっぽい機械学習の勉強もぼちぼち再開します。

この本すごい評判が良いそうで、読んでて構成がわかりやすかったです。 全体構成と今回の対象は下図にします。

f:id:nogawanogawa:20180715133805j:plain:w500

ということで、まずは導入をやってみます。

時系列分析とは

時系列分析とは、時間方向に関して得られる系列的なデータに対する分析のことです。

もうちょっと詳しく言うとこんな感じらしいです。

基本的には関心のある事象における過去・現在・未来の値を適切に把握し(推定し)、関連してその結果を元に、事象の仕組みや影響に関する知見を得たり対策を考えたりする営みである

やることのなんとなくのイメージがこんな感じです。

f:id:nogawanogawa:20180716101546j:plain:w500

時系列分析を通じて、何らかの値を推定することになるのですが、この本では確定的手法と確率的手法の2種類に分けて解説されています。 確定的手法は状態空間モデルを用いた分析(基本)、統計的手法は状態空間モデルを用いた分析(応用)で扱うようです。

確率・統計について

この本とんでもなくわかりやすいので、詳しい説明は本を読んでいただければと。 用語の定義とかをさわりだけ列挙していきます。

平均やら分散やらは、過去にやっているのでこちらをご参照。

tsunotsuno.hatenablog.com

その他、新しく出てきたとこだけ見ていきます。

複数の確率変数の関係

不安定に揺らぐ変数を確率変数と呼び、大文字のアルファベットで表します。 確率変数自体は値ではなく関数に近い意味合いを持ち、確率変数の具体的な値は実現値と言い、小文字のアルファベットで表します。

確率変数X, Yについて、X, Yが同時に成立する同時確率


p(x, y)

と表されます。 また、確率変数X, Yについて、Yが確定したときのX条件付き確率


p(x | y)

と表されます。 条件付き確率に対して、確率の乗法定理を使用すると以下の関係が導けます。


p(x | y)p(y) = p(x, y) = p(y, x) = p(y | x)p(x)

上式から、


p(x | y) = \frac{p(y | x)p(x)}{p(y)}

が得られます。上の関係をベイズの定理と呼びます。 ちなみに、p(x)事前確率p(x | y)事後確率p(y | x)を(xの条件での)尤度(ゆうど)と呼んだりします。事前確率を条件を確定することで事後確率に変換することをベイズ更新と言ったりします。

確率過程

確率変数が時間方向につながった集合を確立過程といいます。 表記としてはX_t, Y_tのように下付き文字を使用してタイムステップを表現します。

共分散・相関

ある時系列を説明する際に使用される統計量として、共分散・相関が使われます。

共分散Cov)は、確率変数のX, Yの関連を表していて、それぞれの確率変数の期待値と比較したときの値の大小の関係性を表現します。


Cov[X, Y] = E[(X - E[X])(Y-E[Y])]

共分散を規格化したものを相関係数\rho)といいます。


\rho = \frac{Cov[X, Y]}{\sqrt{Var[X]}\sqrt{Var[Y]}}

相関係数は-1〜1の値になり、負、0、正に応じて下記のような関係性があることを示しています。

f:id:nogawanogawa:20180716151100j:plain:w500

確率過程X_t, Y_tに関して、時間軸方向にラグkだけずれたX_t, Y_{t-k}に関する共分散を表す相関関数は、下記のように表現されます。


R(t, k) = E[X_t Y_{t-k}]

異なる確率変数間だけでなく、一つの確率変数の時間方向のラグに関する相関を見る場合は、自己共分散・自己相関係数を使用します。

自己共分散は下記の式で表されます。


Cov[X_t, X_{t-k}] = E[(X_t - E[X_t]) (X_{t-k} - E[X_{t-k}])]

自己相関係数は下記の式で表されます。


\rho_{t,k} = \frac{Cov[X_t, X_{t-k}]}{\sqrt{Var[X_t]}\sqrt{Var[X_{t-k}]}}

自己相関係数も-1〜1の範囲で値を取りますが、今回は値ではなくtkの関数になっています。 時刻とラグの関数から周期的な変動を確認することができます。

定常過程と非定常過程

強定常・弱定常・非定常に関して、時刻が変動したときに平均・分散・自己共分散・自己相関係数・確率分布そのものが変動するかどうかをまとめると下記の様になるそうです。

分類 平均・分散・自己共分散・自己相関係数(○:変化する、☓:変化しない) 確率分布(○:変化する、☓:変化しない)
強定常
弱定常
非定常

この本では、弱定常のことを定常と呼んでいるそうです。 要するに、周期性があれば定常、なければ非定常といったところでしょう。

最尤推定ベイズ推定

確率分布を特徴づけるパラメータ\thetaは一般に未知であり、何らかの方法で特定化します。 確率過程Y_t (t=1, 2, ..., T)全体の尤度はp(y_1, y_2, ..., y_T; \theta)で表されます。 この自然対数である対数尤度


l(\theta) = log p(y_1, y_2, ..., y_T; \theta)

で表され、これを用いて尤度が最大にするように\thetaを決定する方法を最尤法と呼びます。 最尤推定


\hat{\theta}= argmax\ l(\theta)

で表されます。

パラメータ\thetaすら確率変数として取り扱うのがベイズ推定と言うそうです。

とりあえずやってみる

数式ばっかりなのは苦手なので、実際にやってみて感じを掴んでみます。 基本的な分析の流れはこんな感じらしいです。

  1. 目的の確認とデータの収集
  2. データの下調べ
  3. モデルの定義
  4. パラメータ値の特定
  5. フィルタリング・予測・平滑化の実行
  6. 結果の確認と吟味
  7. 1.へ戻る

1. 目的の確認とデータの収集

データ

今回使用するデータはこちらです。

  1. ナイル川の年間流量
  2. 待機中の二酸化炭素濃度
  3. 英国の四半期ごとのガス消費量

※本当は4つ目もあったんですが11章に詳しい説明があるらしく、今回は省略。

目的

分析の目的と対象となるデータを以下に示します。

  目的 データ
1 データ取得中に各時点でのノイズをできる限り除去する 1. 2. 3.
2 過去の急激な変化を捉える 1.
3 未来の値を予測する 2.

2. データの下調べ

まずはデータをそのままの形で見てみます。 データの眺め方にもいろいろあるみたいで、この本では下のような表示を行っています。

横軸時間のプロット

ナイル川の流量については、あまり規則性が見られません。 これは不規則な擾乱を除けば、毎年概ね同じ値に保たれているということを意味します。 二酸化炭素濃度については、小刻みに上下を繰り返しながら、全体としては右肩上がりの傾向が見て取れます。 イギリスのガスの量は右肩上がりの傾向があるようで、更に年々変動の幅が大きくなっていることがわかります。

ヒストグラムと五数要約

どのグラフにも大きな外れ値はなさそうです。 ここで、二酸化炭素濃度のデータについては、東日本大震災時の際のデータの欠損がありますので、前後の月のデータの値の平均で代用しています。 欠損は放って置くと変な結果になったりするそうなので、潰しておきます。

自己相関係数

全体として自己相関係数の値がある程度ありそうなので、時系列的な関係性がありそうだとわかります。 また、イギリスのガスのデータはラグが整数のとき自己相関係数が大きくなる傾向があるようです。

周波数スペクトル

フーリエ変換して周波数成分の大きさを確認します。 音とかで言うところの、「聞こえてきた音の中から、成分("ド"とか"レ"とか)がどれくらい入っているか」を調べる作業みたいなものです。

ナイル川の流量についてはあまり周波数成分に特徴が見られません。 二酸化炭素濃度とイギリスのガスについては、特定の周波数成分が大きくなっている事がわかります。

3. モデルの定義

だいたい、傾向がわかったところで仮置きでモデルを設定してみます。

f:id:nogawanogawa:20180717213549j:plain

左から、レベル、傾き、周期を表していて、これを組み合わせでグラフが構成されていそうだということが想像つきます。

  1. レベル
  2. レベル + 傾き + 周期
  3. レベル + 傾き + 周期

なんとなく、こんな感じの組み合わせで表現できそうです。

4. パラメータ値の特定

検討したモデルに対して、具体的な数字を当てはめていきます。 上で指定したモデルでは、経済分野などで使用されるホルト・ウィンタース法が使用できます。

ホルト・ウィンタース法

k期予測値\hat{y}_{t+k}は次のように表されます。


\hat{y}_{t+k}= level_t + trend_t k + season_z{t-p+k_p^+}

このとき、レベル(level)、傾き(trend)、周期(season)はそれぞれ、

level_t = \alpha(y_t - season_{t-1}) + (1 + \alpha)(level_{t-1} + trend_{t-1})

 trend_t = \beta(level_t + level_{t-1}) + (1-\beta)trend_{t-1}

season_t = \gamma(y_t - level_t) + (1 - \gamma)season_{t-p}

となります。ここで、k_p^+ =  \lfloor(k-1)mod p\rfloor +1です。

式で書くとややこしいですが、Rにはこれらをビシッと計算してくれる関数が用意されており、HoltWinters()でパラメータを決定できます。

5. フィルタリング・予測・平滑化の実行

HoltWinters()でサクッとパラメータを特定できます。

赤い点線がモデルによるフィッティング(フィルタリング値)を表すグラフです。

成分ごとに見ていくとこんな感じだそうです。

未来の予測に関しては、二酸化炭素濃度を使用して、未使用の2015年分を予測するとこんな感じになります。

6. 結果の確認と吟味

最後に仮定したモデルが良いか悪いかを評価します。 今回、目的はこんな感じでした。

  • データ取得中に各時点でのノイズをできる限り除去する
  • 過去の急激な変化を捉える
  • 未来の値を予測する

ノイズの除去に関しては、残差(Residual)を確認します。これが大きいとモデルと実測値がかけ離れていることになるので、ノイズを除去しきれていないことを意味します。これはどこまでを許容するかわかりませんが、今回は一旦良しとします。

また、過去の急激な変な変化を捉えるという点については、ナイル川の流量が1899年に急減しているそうですが、そこを捉えることはできていません。この点については更に改良が必要と言うことがわかります。

未来の予測については、二酸化炭素濃度の予測は直感的には一致しているように見えるので、問題ないこととします。

ここで出てきた問題を踏まえて再度手順をやり直し、モデルを改良していきます。

感想

基本的な用語、式の説明と簡単な時系列分析のやり方の説明でした。 久しぶりにやるとRむずいっす。。。