2015年1月31日土曜日

Rでカイ二乗分布の確率密度関数のグラフをシミュレーション的に描く

カイ二乗分布といえば、独立性の検定や適合度の検定で使ったりしますよね。

ウィキペディアには、確率密度関数のグラフが載っていて↓こんな感じ。

カイ二乗分布の確率密度関数のグラフ

このグラフをRで描いてみようと、確率密度関数の式を見てみると、



↑こんなのが載っていて、なんだかすごくむつかしい。Γ(ガンマ)関数ってのが出てきて、さらにそれは積分の形で定義されていたりして。

でもまあ、そのへんは理解できなくても、Rにはカイ二乗分布の確率密度関数(dchisq)が用意されているので、下記のようにすれば、ウィキペディアに載っていたのとそっくりのグラフが描けます。

curve(dchisq(x,1), xlim=c(0,8), ylim=c(0,1), col="black", ylab="dchisq(x, k)")
curve(dchisq(x,2), xlim=c(0,8), ylim=c(0,1), col="blue"   , add=T)
curve(dchisq(x,3), xlim=c(0,8), ylim=c(0,1), col="green"  , add=T)
curve(dchisq(x,4), xlim=c(0,8), ylim=c(0,1), col="red"    , add=T)
curve(dchisq(x,5), xlim=c(0,8), ylim=c(0,1), col="magenta", add=T)
 
#凡例
legend( "topright", lty=1,
        legend = c("k=1"  , "k=2" , "k=3"  , "k=4", "k=5"),
        col    = c("black", "blue", "green", "red", "magenta") )


Rのdchisqを使って描いたカイ二乗分布の確率密度関数のグラフ

で、ここから先は、関数が用意されているにもかかわらず、シミュレーション的なやり方を用いてわざわざ自力で、カイ二乗分布の密度関数曲線を描いてみるという、まったく実用性のない話です。でも、実感みたいなものは湧くかも。

カイ二乗分布の定義はシンプルで、ウィキペディアによると、

独立に標準正規分布に従う k 個の確率変数 X1, ..., Xk をとる。 このとき、統計量

の従う分布のことを自由度 k のカイ二乗分布と呼ぶ。
となっています。

正規分布に従う変数を二乗していくつか足せばいいのか、意外と簡単。

試しに、乱数を発生させて描いてみよう。

n <- 1000000
Z <- (rnorm(n))^2 + (rnorm(n))^2 + (rnorm(n))^2
hist(Z, breaks=seq(0, max(Z)+1, 0.2), freq=F, ylim=c(0,1))




3つ足したんで、k=3の緑の線に対応します。形を見てみると意図通りになってるっぽいですね。

複数のヒストグラムだと重ねたときに見づらいので、カーネル密度推定(density関数)を使って曲線で描画してみましょう。

n <- 1000000
 
X1 <- rnorm(n)
X2 <- rnorm(n)
X3 <- rnorm(n)
X4 <- rnorm(n)
X5 <- rnorm(n)
 
Z1 <- X1^2
Z2 <- X1^2 + X2^2
Z3 <- X1^2 + X2^2 + X3^2
Z4 <- X1^2 + X2^2 + X3^2 + X4^2
Z5 <- X1^2 + X2^2 + X3^2 + X4^2 + X5^2
 
plot(0, 0, type="n", xlim=c(0,8), ylim=c(0,1), xlab="", ylab="")
lines(density(Z1), col="black")
lines(density(Z2), col="blue")
lines(density(Z3), col="green")
lines(density(Z4), col="red")
lines(density(Z5), col="magenta")




k=1、k=2のゼロに近いところでは、かなりずれが大きいですが、それ以外ではけっこう近い曲線になったんじゃないでしょうか。




2015年1月2日金曜日

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.033 μSv/h 0.059 μSv/h  ‐  NA
4 2014/12/31 23:20 0.033 μSv/h 0.059 μSv/h  ‐  NA
5 2014/12/31 23:10 0.033 μSv/h 0.059 μSv/h  ‐  NA
6 2014/12/31 23:00 0.033 μSv/h 0.059 μSv/h  ‐  NA


V5が測定日時時刻、V6が線量の測定値です。(これだけじゃ分かりませんが、V8は高さ1mで換算した推計値らしいです)

V5を横軸、V6を縦軸にすれば、線量の時系列変化が確認できますね。

ただそのままではV5は文字列データなので、日付時刻型に変換する必要がありますね。詳細は以前の記事↓に書きましたが、

Rで文字列を日付時刻型に変換する - Rプログラミングの小ネタ


このデータの場合は as.POSIXlt関数1つで簡単に変換できます。

plot(as.POSIXct(d$V5), d$V8, type="l")


放射線量測定データを時系列プロット(1年分)

さすがに10分ごとに1年分だとデータが多いので、ひと月ぶんだけ取り出してプロットしてみます。

n <- 6 * 24 * 31 # 1ヶ月ぶんの10分の個数
plot(as.POSIXct(d$V5[1:n]), d$V8[1:n], type="l")


放射線量測定データを時系列プロット(ひと月分)

上記は12月ぶんのプロットになっていますね。

細かい変動には意味がないと思われますので、移動平均を取ってみましょう。やり方は過去記事↓にも書きました。

Rで移動平均を求める - Rプログラミングの小ネタ


# 何度か使うので関数化しておく
moving_average <- function(x, n){
  filter(x, rep(1,n)) / n
}
 
時刻   <- as.POSIXct(d_part$V5[1:n])
測定値 <- d_part$V8[1:n]
plot(時刻, moving_average(測定値,3), type="l")


10分前、現時刻、10分後で移動平均を取ったもの

多少見やすくなりましたね。移動平均を取っていないものと、移動平均を取る期間の幅を変えたものを1つのグラフに重ねてプロットし、比較してみます↓

plot( 時刻,                測定値,      type="l", col="gray")
lines(時刻, moving_average(測定値,  6), type="l", col="green")
lines(時刻, moving_average(測定値,144), type="l", col="blue")
legend( "topright",
        lty=c(1,1,1),
        col=c("gray","green","blue"),
        c("測定データそのまま", "1時間の幅で移動平均", "1日の幅で移動平均") )


移動平均なし(グレー)、1時間幅(緑)、1日幅(青)

1日幅の移動平均(青の線)だと、だいぶ傾向が見やすくなったんじゃないでしょうか。山の現れ方には周期があるように見えなくもないです。もし本当に周期性があるなら、それが現れる原因は何か、などを探ってみると面白いかもしれません。




2015年1月1日木曜日

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プログラミングの小ネタ