草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

カーネル密度推定からその後:ksパッケージ

背景

正規分布を逸脱したデータの累積分布曲線を描きたい。

たとえば,こんなデータってことにしよう。

 

f:id:kusanagik:20150612225334p:plain

 

このデータは,

 c(rnorm(100.-5,1),rnorm(100,5,1))->dat

 こうやってつくろう。

Nが200もあれば,quantile関数で,

f:id:kusanagik:20150612225651p:plain

 

で,こんな感じにできる。

でもこれ,飛び飛びだからね。

もし,Nが30程度だったら,

 

f:id:kusanagik:20150612225839p:plain

みたいになっちゃう。

ここで,なめらかな(連続的)曲線,累積分布曲線がほしいって話。

 

経験累積分布関数見てみる

Rには,ecdfっていう関数あるね。

empirical cumulative distribution functionってことね。

まずはこれやっちゃう。えい。

 

ecdf(dat)->ecdfdat

これは,

class(ecdfdat)
[1] "ecdf" "stepfun" "function" 

 っていうことで,関数扱いになる。

なので,

ecdfdat(2)
[1] 0.4074074 

 っていう感じに数字入れてやると累積確率を返してくれる。便利便利。

でも…

plot(ecdfdat)

f:id:kusanagik:20150612230608p:plain

っていう感じになっちゃうの。Nが小さいと。

さて。

 

とりあえずカーネル密度推定

ま,とりあえずカーネル密度推定で確率密度曲線を得よう。

カーネル密度推定の関数っていっぱいあるよね。

densityでもいいけど,ここではksパッケージ使おう。

ksパッケージいいよ,好きよ。

 

http://cran.r-project.org/web/packages/ks/ks.pdf

 

library(ks)

 

カーネル密度推定はバンド幅とかが命なのだけど,ちょっと端折って,えい。

これだけ。

 kde(dat)->kdedat

plot(kdedat)

 f:id:kusanagik:20150612231651p:plain

まあいいとしよう。

積分布曲線描きたい

さて,ここでこのカーネル密度推定による確率密度関数を,数値積分繰り返していって累積分布曲線描くのか…って考えた。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")

f:id:kusanagik:20150612233229p:plain

 

お,よさげよさげ。

カーネル密度推定で得た確率密度関数に従う乱数発生

これすごい。rkdeっていう関数。

rkde(fhat=kdedat,n=100)

nに乱数の数指定してやる。

ためしに1,000くらい乱数作ってまたkdeしてプロットしちゃえ。

f:id:kusanagik:20150612233832p:plain

おお,いいねいいね。平滑化ってかんじよね。あと,これで結構シミューレーションとか自由になるね。

この乱数を使う?

パラメトリックブートストラップみたいな感じで,中央値の信頼区間とかやってみる? 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  

こんな感じか。

まあ,でもバンド幅次第だよなあ…