カーネル密度推定を君へ

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