論理の流刑地

流罪に遭い戸惑う世捨人の雑記

カーネル密度推定を君へ

Motivation

数値ベクトルxが与えられているとき
xの値が与えられていない区間についても、なめらかな形で確率密度関数を求めたいときがある。

そんなときに多く用いられるのが、カーネル密度推定である。
カーネル関数をK(y)を用いた以下の式で、(実際の観測値の有無にかかわらず)任意のf(X = x)の値を求めることができる。

 f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n K(\frac{x - x_i}{h})

しかし、Rのstats::densityは実はこの式を直接実行することはできない。
(from,to,n などの引数をいじって無理やり推定点をずらしても、内部ではapproxを用いて近傍の2つの点の補間を行なっているので、上記の式が直接適用されているわけではない)
ksパッケージのkdeも同様で、内部でapprox()を使っているので直接上記の式を計算しているわけではない。
(結構な落とし穴だと思うがわかりやすく明言されていない)

なので、カーネル密度推定を用いた統計手法を用いる時は、自分で上記式の計算関数を実装するのがよい。

ガウス関数の場合

カーネル関数としてガウス関数(標準正規密度関数)を用いるとすると、
 f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n K(\frac{x - x_i}{h})
 f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n \frac{1}{\sqrt{2\pi}}exp({-\frac{(\frac{x-x_i}{h})^2}{2}})
 f(x) = \frac{1}{nh}\Sigma_{i = 1 }^n h \times \frac{1}{\sqrt{2\pi h^2}}exp({-\frac{({x-x_i})^2}{2h^2}})
 f(x) = \frac{1}{n}\Sigma_{i = 1 }^n \frac{1}{\sqrt{2\pi h^2}}exp({-\frac{({x-x_i})^2}{2h^2}})

ここで\Sigmaの中を見てみると、これは平均x_i、分散がh^2正規分布の密度関数となっている。
だから、所定の方法でband幅hさえ求めてしまえば、あとはdnorm(x,x_i, h)の平均を求めればそれが密度となる。
各ケースにweightが付与されている場合は、平均が重み付き平均になるだけだ。

海外サイトのMLでも実装例が示されている(weightには対応してない)

実装

Rで実装する。weightにも対応するようにしている。
バンド幅は、bw.SJ()とbw.nrd0()関数を利用して求めている

density_pntEst <- function(x,  at , w=NULL ,bw = "sj-ste"){
# list-wise deletion
if(is.null(w)){
w <- rep(1 , length(x))
}
df <-  data.frame(x=x , w = w)
idx <- complete.cases(df)
x <- df[ complete.cases(df) , "x"]
w <- df[ complete.cases(df) , "w"] 
bw <- tolower(bw) #小文字化

# band width
if(bw %in% c("sj", "sj-ste")){
bwval <- bw.SJ( x , method = "ste")
}else if( bw == "sj-dpi"){
bwval <- bw.SJ( x , method = "dpi")
}else if(bw == "nrd0"){
bwval <- bw.nrd0(x)
}else{
stop("band width決定手法が適切ではない")
} #if

w <-  (w * length(x) )/ sum(w) # weightを平均1に調整

at <- matrix(at , ncol = 1) 
y <- apply(at , 1, function(a,x ,bw,w ){
return( mean(  dnorm( a, x , bw)  * w )) 
} , x= x , bw = bwval,w =w)

return( list(x = x , y = y , bw = bwval , index =idx  ) )

} #function

汚いコードなのは重々承知だが、何百回も呼び出すわけじゃないからご愛嬌。

www.youtube.com