二次元カーネル密度推定

カーネル密度推定とは標本値をぼやけさせて確率密度関数を連続的に推定することで最頻値を特定する手法のひとつである。"kde2d"を用いた2次元カーネル密度推定を行う。"2次元"と表記されているものの、扱うデータは2次元でも3次元でも良い。

以下のようなある地域の一年間における犯罪が起こった26地点の緯度および経度が記されているデータがあるとする。このデータからどの地点、どの辺りが最も犯罪発生率が高く危険な地点かをカーネル密度推定にて特定する。

kde2dデータ

まず、以下のようにエクセルから緯度カラムのデータをコピーする。

kde2dデータのコピー

コピーをしたら、そのままRを起動させ(または起動させて置いたRで)、以下のコマンドを打ち、"x"に緯度データを取り込む。

1|$x=scan("clipboard")

同様に、経度カラムのデータをコピーし、"y"に取り込む。

1|$y=scan("clipboard")

まずはXY平面にプロットしてみる。この操作はカーネル密度推定には直接は関係ない。

1|$plot(x,y,xlab="latitude",ylab="longitude")
プロット

点が集中しているところが犯罪発生率が高いところであるが、カーネル密度推定を行うとピンポイントで最も犯罪発生率が高い座標を特定できる。カーネル密度推定を実行するためのコマンド"kde2d"を使用するためにパッケージ"MASS"をインストールして読み込む。既にインストールしてる場合は、1行目のコマンドは不要。

1|$install.packages("MASS", repos="http://cran.ism.ac.jp/")
2|$library(MASS)

次に、カーネル密度推定を行う。"kde2d"の使い方は以下のようになっている。

kde2d(x, y, h, n, lims = c(range(x), range(y)))
"x"は緯度データ(上で作成した変数"x")
"y"は経度データ(上で作成した変数"y")
"h"はカーネル密度推定をする際の各々の軸に対するバンド幅を決定するベクトルであり、"i"および”j"を数値とすると"c(i,j)"の形で指定する。バンド幅は、自身のデータに合わせて一番良い結果が得られるように適当な値を見つけ出す。これを補助してくれるのがページ下部にまとめたバンド幅決定関数である。今回は"bandwidth.nrd"を使用する。
"n"はカーネル密度推定をする際に、どれだけのグリッドを発生させるか、を示す。この値を大きくすると結果として得られる座標がよりきめ細かになる。
"lims"は"xi", "xj", "yi", "yj"を数値とするとc(xi, xj, yi, yj)のように入力し、各座標軸で密度推定をする範囲を指定する。例えば、c(-12,0,-14,4)のようにするとx軸の密度推定は-12から0、y軸の密度推定は-14から0の範囲でのみ行われる。基本的には指定しなくて良い。

以上を踏まえて、以下のコマンドで"kde2d"を実行し"d"に結果を格納する。

1|$d=kde2d(x,y,c(bandwidth.nrd(x),bandwidth.nrd(y)),n=80)

以上の計算で"d"に主に以下のデータが格納される。
"d$x"に密度計算後のx座標
"d$y"に密度計算後のx座標
"d$z"に対応する(x,y)の密度

得られた結果を以下のコマンドで画像化する。白く光っているところが犯罪発生件数が多い場所である。

1|$image(d,xlab="latitude",ylab="longitude")
イメージ

以下のコマンドで等高線を描く。

1|$contour(d,xlab="latitude",ylab="longitude")
等高線

最後に、"d$z"の最大値を以下のコマンドで検索する。

1|$max(d$z)

得られた数値 (0.005004014)を"d$z"から検索し、そのときの(x,y)が最高密度座標である。結果として(x, y)=(4.759494, 4.544304)が最高密度座標であった。このことから、緯度4.759494および経度4.544304の地点が最も犯罪発生率が大きい地点であると特定された。等高線の頂点も大体その位置を示していることがわかる。

-----バンド幅を決定する関数には以下のものがある。

FunctionPackage
bandwidth.nrdMASS
ucvMASS
bcvMASS
width.SJMASS
bw.nrd0stats
bw.nrdstats
bw.ucvstats
bw.bcvstats
bw.SJstats

-----3次元データの扱いについて

また、このようなデータ形式ではなく、以下のようなカラム1に緯度、カラム2に経度、カラム3に犯罪発生件数が示されている場合、つまり3次元データもkde2dで扱える。

3次元データ

この場合は、以下のように、緯度データおよび経度データを対応する犯罪発生件数分複製する。(緯度, 経度)=(2, 5)のときの犯罪発生件数が6なら(2, 5)を6行作成する。これを全ての座標で行う。ちょっとしたプログラムを書くと早い。あとは、緯度データを"x"に、経度データを"y"に格納し、上と同じ計算である。これは各座標をグリッドとみなし、犯罪発生件数のデータがないところは、犯罪発生件数が0として、緯度および経度の座標を0個発生させているとみなせる。

3次元データ入力

このようにkde2dを用いることで、3次元データを扱えるようになるが、これは"kde3d"にもあてはまり、"kde3d"では(x, y, z)座標とその座標が持つ何らかの値からなる4次元のデータを扱えるようになる。

このエントリーをはてなブックマークに追加

Site search

ページのトップへ戻る