Rで移動平均を求める

Rには移動平均そのものずばりを求める関数はないようです。でも、filter関数を使えば簡単に実現できます。

まずは、filter関数の動きを確認。

xを時系列データとします。

> x <- c(1, 2, 3, 1, 2, 9, 4, 2, 6, 1)
> x
 [1] 1 2 3 1 2 9 4 2 6 1

このxにfilter関数を適用すると、

> filter(x, c(1,1,1))
Time Series:
Start = 1
End = 10
Frequency = 1
 [1] NA  6  6  6 12 15 15 12  9 NA

となります。

c(1,1,1)の3つの要素は、（時点 n-1，現時点 n，時点 n+1）に対応します。重みに差を付けられるのですが、ここでは同等にしたいので、すべて1にしています。端の要素は前後の要素がないのでNAとなります。

↓こんな感じですね。

Rのfilter関数の適用イメージ

これだけだと、（時点 n-1，現時点 n，時点 n+1）を足したものになってしまうので、足した要素の数（上記だと3）で割れば、移動平均を取ることになりますね。

ということで、ごくシンプルな移動平均は↓こんな感じで実現できます。

filter(x, c(1,1,1)) / 3

平均を取る幅を広げて、n-2, n-1, n, n+1, n+2 のように5つ分にしたい場合は、

filter(x, c(1,1,1,1,1)) / 5

とやればOK。

分かりやすいように c(1,1,1)、c(1,1,1,1,1))とベタに書きましたが、数が増えた場合も考えて rep(1,3)、rep(1,5)のように書くのが普通ですかね。

何度も使うなら、

moving_average <- function(x, n){
  filter(x, rep(1,n)) / n
}

のように関数化してしまってもいいかもしれません。

【蛇足】
最初「filter」っていう関数名がピンとこなかったんですが、画像処理で使うフィルタと同じ概念なんだと気付いて納得しました。あるピクセルの周囲のピクセルと平均を取れば、画像は少しぼんやりするけど、細かいノイズを消したり軽減したりすることができるという、あの類の処理です。移動平均ってローパスフィルタなんですね。

[2015年1月2日追記]
実際のデータの時系列プロットで使ってみたサンプルです↓
Rで放射線モニタリング情報を時系列プロットする - Rプログラミングの小ネタ
■