Rプログラミングの小ネタ

データを公開するときはCSV形式みたいに、RAWデータかそれに近い形で提供してくれると、いろいろと利用できて嬉しいですよね。 原子力規制委員会のサイト では、モニタリングポストで測定された空間線量のデータをCSV形式で公開しているので、これを使って時系列データの可視化を試してみましょう。 下記のページで、都道府県やエリア、期間の日付を選択してダウンロードできます。 ダウンロード | 東日本大震災関連情報 放射線モニタリング測定結果等 | 原子力規制委員会 試しに東京都の都健康安全研究センターのモニタリングポストのデータを2014年1月1日～2014年12月31日の1年分ダウンロードしてみました。10分ごとのデータなのでかなりの量ですね。約5.5MBくらいありました。 まずは、読み込んで、冒頭部分を表示。 > d <- read.csv ( "新宿2014.csv" , header=F , as.is=T ) > head ( d ) V1 V2 V3 V4 1 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 2 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 3 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 4 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 5 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 6 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 V5 V6 V7 V8 V9 V10 V11 1 2014 / 12 / 31 23 : 50 0.034 μSv/h 0.061 μSv/h ‐ NA 2 2014 / 12 / 31 23 : 40 0.033 μSv/h 0.059 μSv/h ‐ NA 3 2014 / 12 / 31 23 : 30 0....

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」っていう関数名がピンとこなかったんですが、画像処理で使うフィルタと同じ概念なんだと気付いて納得しました。ある...