カーネル密度推定

密度を推定する

　データが点の集まりとして手に入ったとしましょう(例えば身長・体重データとか喫煙歴とか)。それで、この点々の母集団の分布を考えるときは、例えばこの点々を基にしてQQプロット打ち込んでう～ん正規分布っぽいから平均と分散を推定するか！　とかやっていくわけです。これがパラメトリックな密度の推定になります。
　さて、ノンパラメトリックな密度の推定を考えてみましょう。つまり、母集団の分布が属する"族"を(ほとんど)設定せず、ただ単に今手元にある点々から、それっぽい密度関数をつくりたい……そういう需要はあると思います。めっちゃあります。
　というわけで、今回はノンパラメトリックな密度推定法の、カーネル密度推定と呼ばれる手法について軽くまとめ&実装してみました。

　1次元の場合を考えます。
　カーネル密度推定において、点x上での密度の推定値は、バンド幅をh、サンプルサイズをnとしたときに、

$\displaystyle{ \hat f (x) = \frac{1}{nh} \sum _{i=1}^{n} K\left( \frac{x-X_i}{h}\right) }$

によって表されます。ここで、 $\displaystyle{ K() }$ はカーネル関数で、じゃあカーネル関数ってなんだよというと、

$\displaystyle{ \int_{-\infty}^{\infty}K(t)dt = 1 }$

なる関数のことです。
　バンド幅ってなんだよとなると思いますが、まあイメージとしてはヒストグラムのビン幅みたいなものです。このパラメータhのことをバンド幅って呼ぶんだなあぐらいの認識でも良いと思います。
　なお、通常、このカーネル関数は原点対称なものが使われ、例えば正規分布の密度関数なんかはガウシアンカーネルと呼ばれています。
　このカーネル関数を色々と設定することによって、ただのサンプル(点の集まり)から密度推定が可能になるのです。

　カーネル関数の定義より、この推定された密度関数がきちんとマイナス無限大から無限大で積分したときにその値が1になっているとわかりますね。

※統計学でカーネルというと文脈によって複数意味があるので注意。ここでいうカーネル関数というのはノンパラメトリック統計におけるカーネル関数のことです。ベイズ統計や機械学習を勉強したことある人なら実感したことあると思います。参考→カーネル(統計学)　ちなみにマルコフ過程の推移核のことをマルコフカーネルとか確率カーネルとか言ったりもする(M(x,dy)みたいな表記)。　個人的に、カーネルって言われても何のこと言ってるのかよく分からないので正直この言葉使ってほしくないですね。

　カーネル関数の種類は多分上げてったらキリがないと思うんですけど、さっきのWikipediaの記事にいくつかまとまっているので気になって人は読んでみてください。

非心t分布の密度関数を図示したい

$\displaystyle{ U \sim N(\lambda , 1) \\ V \sim \chi ^2(m) \\ UとVは独立 }$

という状況を考えましょう。このとき、

$\displaystyle{ \Large \frac{U}{\sqrt{\frac{V}{m}}} }$

が従う分布のことを、非心度λ、自由度mの、非心t分布と呼びます。t(m,λ)で表すことにしましょうか。例えば仮説検定で母集団同士の平均を比べるときに平均に差があるような場合はt統計量がこの分布に従いますね。

　この定義から明らかに、非心度λが0である非心t分布、t(m,0)は、自由度mのt分布と一致することがわかります。
　さて、ではt(m,λ)の密度関数はどんな関数になっているかというと、以下のような関数になっています。(求め方は普通にt分布の密度求める時と一緒なので、見た目ほど難しくはないのですが、計算はめんどくさいです)

$\displaystyle{ \LARGE f(t\ ;m,\lambda)=\frac{exp(-\frac{\lambda ^2}{2})}{\sqrt{\pi m} \ \Gamma(\frac{m}{2}) } \sum _{j=0}^{\infty} \frac{\lambda ^j t^j \Gamma \left(\frac{m+j+1}{2}\right)2^\frac{j}{2}}{j\ ! \ m^{\frac{j}{2}} \ \left(1+\frac{t^2}{m}\right)^{\frac{m+j+1}{2}}} }$

なんだこれ！？

　さて、私はちょっと気になったんですけど、t(m,λ)って、λの値を動かしたときにどういった感じで形が変わっていくんでしょうか？　定義より明らかにλを動かすほど分布が横にスライドしていきそうなものですが……ほかには、expの係数から、λが0より離れていくほど、分布が潰されていきそうな感じがするでしょうか。まあ、ぱっと見じゃ、そのくらいのことしかわからないですね。

　ということで、非心t分布の密度関数を私は図示したいのです。お気楽に、お手軽に。
　この関数には無限和が含まれているのでどこかのタイミングで計算を打ち切って密度関数を書くしかないですが、もう一つ、密度関数を描く方法がありそうですね。

　この分布は式は複雑ですが、この分布に従う乱数の発生は比較的容易ですから、もう乱数を大量に発生させてしまって、そのサンプルから密度をノンパラメトリック推定してしまおう～！　というのが、目標です。

　じゃあまず乱数を沢山作りましょう。自由度は10にします。非心度は今はとりあえず5にしておきます。後でいじります。

n <- 10000
m <- 10
lambda <- 5
U <- rnorm(n=n , mean = lambda ,sd = 1)
V <- rchisq(n=n , df=m)
RANSU <- U / sqrt(V/m)
hist(RANSU)

f:id:hasegawa-yuta289:20210814021612p:plain

　さて、これによってt(10,5)に従う乱数が10000個できました。
　ヒストグラムだけでも、まあ正直なんとなく密度関数の形が見えてきますけど……
　ここから、カーネル密度推定を、上でやったように行いましょう！

　バンド幅hは色々と設定できるので、とりあえずh=0.2にしておきます。
　カーネル関数にはガウシアンカーネルを利用しましょう。標準正規分布の密度関数なので、以下のような関数になっています。

$\displaystyle{ K(t) = \frac{1}{\sqrt{2\pi}}exp \left(-\frac{t^2}{2}\right) }$

再掲：カーネル密度推定

$\displaystyle{ \hat f (x) = \frac{1}{nh} \sum _{i=1}^{n} K\left( \frac{x-X_i}{h}\right) }$

これを参考にしてプログラミングしていきましょ～

h <- 0.2
f.hat.Gaussian <- function(x,h,X){
  n <- length(X)
  ATAI <- c()
  for(i in 1:length(x)){
  ATAI[i] <- sum(dnorm((x[i]-X)/h)) / (n*h)
  }
  return(ATAI)
}
hist(RANSU , freq=FALSE , main = paste("ガウシアンカーネルによる密度推定 h=" , h , sep="") , 
     xlab = paste("x (n=" , n , ")", sep=""))
curve(f.hat.Gaussian(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col="red" , lwd=2 , add=T)

f:id:hasegawa-yuta289:20210814021633p:plain

ということで、ガウシアンカーネルを使った非心t分布の密度関数の推定ができました。
　nを増やし、hを細かくしていけばより精緻な推定ができるでしょう。

　ちなみに、試しに三角カーネルを利用したときの密度推定は以下のようになります。Iは定義関数です。(その他の設定は全部同じ)

$\displaystyle{ 三角カーネル:K(t) = (1-|t|)I(t\in [-1 , 1]) }$

h <- 0.2
f.hat.triangular <- function(x,h,X){
  n <- length(X)
  ATAI <- c()
  for(i in 1:length(x)){
    t <- (x[i] - X) / h
    ATAI[i] <- sum( (1 - abs(t)) * (abs(t)<1) ) / (n*h)
  }
  return(ATAI)
}
hist(RANSU , freq=FALSE , main = paste("三角カーネルによる密度推定 h=" , h , sep="") , 
     xlab = paste("x (n=" , n , ")", sep=""))
curve(f.hat.Gaussian(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col="red" , lwd=3 , add=T)
curve(f.hat.triangular(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col="green", lwd=2 , add=T)
legend("topright" , legend=c("ガウシアンカーネル","三角カーネル") , lty=1,col=c("red","green"),lwd=2)

f:id:hasegawa-yuta289:20210814021646p:plain

　まあ、大体同じ結果になりましたね。実用上はどっちを使っても問題なさそうですが、どうでしょう？

　さて、ここから本題です。t(m,λ)のλを色々動かしたときの、非心t分布の変化を見ていくことにします。
　バンド幅は0.05で、n=1000000にしましょう。自由度は10として、非心度を0,2,4,……,20と変化させます。

curve(dt(x,df=10),xlim=c(-2,24) , col="red" , lwd=2 , main="非心t分布(自由度10 , 非心度0～20)",ylab="密度")
n <- 1000000
m <- 10
h <- 0.05
for(lambda in seq(2,20,2)){
  U <- rnorm(n=n , mean = lambda ,sd = 1)
  V <- rchisq(n=n , df=m)
  RANSU <- U / sqrt(V/m)
  curve(f.hat.Gaussian(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col=rgb(red=(20-lambda)/20 , blue=lambda/lambda , green=0) , lwd=2 , add=T)
}

f:id:hasegawa-yuta289:20210814021701p:plain

ということで、t(10,λ)の密度関数は、λを0から20まで動かすとこの図のようになることがわかりました。
　分布が右にスライドすることや、段々と分布が潰れていくことは、予想通りの結果になりましたね。

　ということで、今回はサンプルから密度を推定する、カーネル密度推定を行いました。ちなみに今回作ったRのコードをいちいち書かなくても、Rにはこれと同じことをやってくれるdensity()という関数があるのでそれを使えば簡単です。便利ですね。ちなみにバンド幅hの良い選び方とかもあるんで、気になった人は↓の本を読んだり調べたりしてみてください。(おわり)

参考

Ｒによる計算機統計学

作者:Maria L..Rizzo
オーム社

Amazon