カーネル密度推定からその後:ksパッケージ
背景
たとえば,こんなデータってことにしよう。
このデータは,
c(rnorm(100.-5,1),rnorm(100,5,1))->dat
こうやってつくろう。
Nが200もあれば,quantile関数で,
で,こんな感じにできる。
でもこれ,飛び飛びだからね。
もし,Nが30程度だったら,
みたいになっちゃう。
ここで,なめらかな(連続的)曲線,累積分布曲線がほしいって話。
経験累積分布関数見てみる
Rには,ecdfっていう関数あるね。
empirical cumulative distribution functionってことね。
まずはこれやっちゃう。えい。
ecdf(dat)->ecdfdat
これは,
class(ecdfdat)
[1] "ecdf" "stepfun" "function"
っていうことで,関数扱いになる。
なので,
ecdfdat(2)
[1] 0.4074074
っていう感じに数字入れてやると累積確率を返してくれる。便利便利。
でも…
plot(ecdfdat)
っていう感じになっちゃうの。Nが小さいと。
さて。
とりあえずカーネル密度推定
ま,とりあえずカーネル密度推定で確率密度曲線を得よう。
カーネル密度推定の関数っていっぱいあるよね。
densityでもいいけど,ここではksパッケージ使おう。
ksパッケージいいよ,好きよ。
http://cran.r-project.org/web/packages/ks/ks.pdf
library(ks)
カーネル密度推定はバンド幅とかが命なのだけど,ちょっと端折って,えい。
これだけ。
kde(dat)->kdedat
plot(kdedat)
まあいいとしよう。
累積分布曲線描きたい
さて,ここでこのカーネル密度推定による確率密度関数を,数値積分繰り返していって累積分布曲線描くのか…って考えた。integrate(kdedat, -Inf,.001)を1まで繰り返すかんじっしょ…ってやったけど,kdeが返すオブジェクトは関数じゃないんかいww
ていうかksパッケージには,pkdf(累積分布関数)っていう関数あるんだ。
pkde(fhat=kdedat,q=seq(min(dat),max(dat),.01))
これをプロットしたらいいね。
ran<-seq(min(dat),max(dat),.001)
plot(ran,pkde(fhat=kdedat,q=ran),type="l",xlab="Value",ylab="p")
お,よさげよさげ。
カーネル密度推定で得た確率密度関数に従う乱数発生
これすごい。rkdeっていう関数。
rkde(fhat=kdedat,n=100)
nに乱数の数指定してやる。
ためしに1,000くらい乱数作ってまたkdeしてプロットしちゃえ。
おお,いいねいいね。平滑化ってかんじよね。あと,これで結構シミューレーションとか自由になるね。
この乱数を使う?
パラメトリックブートストラップみたいな感じで,中央値の信頼区間とかやってみる? B = 1,000でパーセンタイル法でいいか。
result<-numeric(0)
for(i in 1:1000){
result[i]<-median(rkde(fhat=kdedat,n=1000))
}
quantile(result,c(.025,.975))2.5% 97.5%
2.733980 3.570161
こんな感じか。
まあ,でもバンド幅次第だよなあ…