Rで"Kernel Density Estimation"、その方法と本質!?
大学では、位置情報関連のデータの分析を主にしてます。
そこでいつもよく使う“カーネル密度推定”なんですが、GUIでサクッとやってしまうとその本質を学ばずして、してやったりしがちなので、ちょっと詳しく参考書を読みながら学び直し！
カーネル密度推定とは...
“単変量または多変量確率密度関数を推定する手法...
観測されたデータを内挿し、平滑化された推定値を得ることができる。”*1
対象地域R上に観測されたイベントの地点sがn個あるとし、
カーネル関数及びバンド幅を決定すれば、カーネル密度は以下の数式で推定できる。
kのカーネル密度関数は、原点に対して対象な2変量確率密度関数、
τはカーネルバンド幅で、この値は常に正の値をとり、平滑化の程度を調整する。
δは、エッジ補正を目的とした係数。（δの数式は割愛）
このカーネル関数kやカーネルバンド幅rを式に与えて、λを細かいグリッド上で繰り返し計算することで、対象地域Rにおける強度λの変動パターンを視覚化できる。
ということだった。
ちょっと小難しい数式のわがんねぇような話は、このくらいにして、
以下からRでのちょこちょこやってみよう。
まずは、τ=1とした場合のカーネル関数（4次関数）の値をプロット！
u <- (0:100)/100 tau <- 1 # 以下の計算は、カーネル関数（四次関数）の式に基づく k <- 3/(pi * tau^2) * (1 - u^2/tau^2)^2 plot(x, k, type = "l", xlab = "distance", ylab = "")
↑はそのx軸の距離が大きくなるに従って、y軸の値が減少していくということは、距離が遠くなるに従って観測地点sの評価は低くなるってことを示してる。
つまり、
バンド幅が大きいと、全体的的傾向は観察できるが局所的変動を観察できず、
バンド幅が小さいと、全体的傾向は観察できないが、局所的変動を観察できる。
とのことだろうか。後日、検証してみよう。
（このカーネル関数の代表例としては、ガウス関数、イパネクニコフ関数、四次関数（biweight関数、quartic関数）が挙げられ、今回は四次関数を利用。）
次にRにプリセットで備わっているbodminデータを使って、可視化までやってみよう。（bodminデータは花岡岩の岩山の地点データ。）
⇒自分のデータも加工したものであれば公開してよいはずなので、今度！！
# カーネル密度推定を実行するpkg library(splancs) # 実験データ data(bodmin) # plot呼び出し plot(bodmin$poly,asp=1,type="n",axes=F,xlab="",ylab="") # エリアデータ polymap(bodmin$poly,add=T) # カーネル密度推定の可視化、hoがカーネルバンド幅だよ。 image(kernel2d(as.points(bodmin),bodmin$poly,h0=2,nx=100,ny=100),add=T,col=terrain.colors(20)) # ポイントと行政区域の追加 pointmap(as.points(bodmin),add=T)
こんなふうに可視化できるみたい。
バンド幅を変えた比較等も行えるので、近いうちにやってみよう。
このバンド幅τによって、視覚的な結果も全然変わってくるということでで、可視化する際は注意が必要やな。
そのバンド幅を調整する方法としては、単純に見た目で決めても良いそうだし（ホンマに！？）、
ちゃんと計算するには、データ数Nを標準偏差または四分位範囲IQRの小さい方を使う方法や、
最小二乗誤差法、交差検証対数尤度関数（cvloglk関数）等で算出する方法が存在しているので、
カーネル密度推定を使うときは、このバンド幅の調整をしっかりとやることが重要ですね！！
可視化だけやとすぐできちゃうので。笑
今回参考にした著書は、
地理空間データ分析 (Rで学ぶデータサイエンス 7)
Ｒによる空間データの統計分析 (統計科学のプラクティス)
です。
次回は、バンド幅の比較と最適な推定方法かマーク付き点過程（多変量ってことかな）、他にはRTB関連の論文のレビューのどれかでもしようかな。
