三次元カーネル密度推定

二次元カーネル密度推定と同様に三次元カーネル密度推定も最頻値を特定する手法のひとつである。"kde3d"を用いた3次元カーネル密度推定を行う。"3次元"と表記されているが、扱うデータは4次元である。この"kde3d"を用いると四次元カーネル密度推定も可能である。もちろんシンプルに3次元データを扱うこともできる。

以下のように、Xは (-14, -12, -10, -8, -6, -4, -2, 0)の8点、Yは (-8, -6, -4, -2, 0, 2, 4)の7点、Zは (-10, -8, -6, -4, -2, 0)の6点によって336点 (=8・7・6)の格子点が構成されており、それぞれの格子点が0以上80以下の整数からなるある値Sを持つとき、Sが最大となるXYZ空間中の一点を特定する。もちろん336点の格子点は実際の値を持っているので最高値は一意に決まるはずであるが、それは336個の格子点の中の最大点である。もし、2刻みの格子点ではなくさらに細かい刻みでXYZ空間中における最高点を推定したいとき、また、336点の実測値における最高点はデータのばらつきによって特異的に生じた最高点である可能性もあるが、そのような値を嫌う場合には3次元カーネル密度推定が有用である。

kde3d用のデータ

まず、以下のグラフにて3次元空間に格子点とその値を色彩強度にて示した。赤い点で示されているところが大きな値を示す点 (全格子点の10%)である。青い点はその次に大きな値 (全格子点の30%)とする。この空間中の最高点を特定する。

データを三次元散布図にしたもの

まずは、以下のようにそれぞれの格子点をSの値分だけ複製する。今回はSは0以上80以下の整数であったのでこのようなことができるが、Sの値が小数を含む場合であったり、小さすぎたり、大きすぎたりする場合はSの値は扱いやすい整数になるように適当に加工する。ただし、データをそのように扱っていいと判断できるときのみである。以下の作業は手間がかかるのでbashとかperlか何かで簡単なスクリプトを書くと良い。

データをSの分だけ複製したもの

以下のようにXの値を選択しコピーする。

コピー選択

そのままRのコンソールに移り、以下のコマンドを打つことで変数"x"に値を格納する。

1|$x=scan("clipboard")

"y"および"z"も同様のコマンドにて作成する。

1|$y=scan("clipboard")
2|$z=scan("clipboard")

以上でインプットデータが揃ったので、以下のコマンドで2つのパッケージを読み込む。パッケージのインストールはこちら。"misc3d"はカーネル密度推定に用いる。また、"MASS"はバンド幅を決定するために用いる。

1|$library(misc3d)
2|$library(MASS)

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

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

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

1|$d=kde3d(x,y,z,c(bandwidth.nrd(x),bandwidth.nrd(y),bandwidth.nrd(z)),n=40)

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

得られた結果を以下のコマンドで画像化する。

1|$image3d(v=d$d,x=d$x,y=d$y,z=d$z,)

白く光っているところが密度が高いところ、すなわちこの場合であると値Sが大きな領域である。

イメージ3D画像

しかし以上では軸もなく分かりにくいのでgnuplotにて画像化した。色彩強度で表されており、紫が密度が低い領域、黄が密度の高い領域を示している。以下はy=-1.5付近にてXZ平面で切り取った断面図である。

4次元イメージ画像

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

1|$max(d$d)

得られた数値 (0.0004675675)を"d$d"から検索し、そのときの(x, y, z)が最も密度が高い点、この場合であると値Sが最も大きい点である。結果として(x, y, z)=(-8.974359, -1.538462, -4.358974)が最高点であると特定できた。これは、gnuplotによる4次元散布図と照らし合わせても確からしいと判断できる。ちなみにx軸、y軸、z軸および密度をファイルに落とすときは以下のように上の"n"に合わせて改行すると見やすい。

1|$write(d$d,"./d.txt",ncolumns=40)
2|$write(d$x,"./x.txt",ncolumns=40)
3|$write(d$y,"./y.txt",ncolumns=40)
4|$write(d$z,"./z.txt",ncolumns=40)

バンド幅を決定する関数については二次元カーネル密度の項に示した。

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

Site search

ページのトップへ戻る