多変量カーネル密度推定

カーネル密度推定は、統計学における基本的な問題の 1 つである密度推定、つまり確率密度関数の推定のための ノンパラメトリック手法です。これは、統計特性が改善されたヒストグラム密度推定の一般化と見ることができます。ヒストグラムの他に、他の種類の密度推定には、パラメトリック、スプライン、ウェーブレット、フーリエ級数などがあります。カーネル密度推定は、1950 年代と 1960 年代に単変量データに関する科学文献で初めて導入され^{[1] [2]}、その後広く採用されました。多変量データに対する類似の推定が多変量統計学への重要な追加となることがすぐに認識されました。1990 年代と 2000 年代に行われた研究に基づいて、多変量カーネル密度推定は、単変量カーネル密度推定に匹敵する成熟度に達しています。^{[3] [4] [5]}

ヒストグラムの構築方法を説明するために、50 点の合成 二変量データセットを例に挙げます。このためには、アンカー ポイント (ヒストグラム グリッドの左下隅) を選択する必要があります。左側のヒストグラムでは (-1.5, -1.5) を選択し、右側のヒストグラムでは、アンカー ポイントを両方向に 0.125 ずつシフトして (-1.625, -1.625) にします。両方のヒストグラムのビン幅は 0.5 なので、違いはアンカー ポイントの変更によるものだけです。色分けは、ビンに含まれるデータ ポイントの数を示します。0 = 白、1 = 淡黄色、2 = 明るい黄色、3 = オレンジ、4 = 赤です。左側のヒストグラムでは、上半分の密度が下半分よりも高いように見えますが、右側のヒストグラムではその逆であり、ヒストグラムがアンカー ポイントの配置に非常に敏感であることが確認できます。^[6]

このアンカー ポイント配置の問題に対する 1 つの解決策は、ヒストグラムのビニング グリッドを完全に削除することです。下の図の左側では、カーネル (灰色の線で表示) が上記の 50 個のデータ ポイントのそれぞれを中心としています。これらのカーネルを合計した結果が右側の図に示されており、これがカーネル密度推定値です。カーネル密度推定値とヒストグラムの最も顕著な違いは、前者にはビニング グリッドによって誘発される人工物が含まれていないため、解釈が容易なことです。色付きの等高線は、それぞれの確率質量を含む最小の領域に対応しています。赤 = 25%、オレンジ + 赤 = 50%、黄色 + オレンジ + 赤 = 75% であり、単一の中央領域に最も高い密度が含まれていることを示しています。

密度推定の目的は、有限のデータ サンプルを取得し、データが観測されていない場所も含め、あらゆる場所で基礎となる確率密度関数について推論することです。カーネル密度推定では、各データ ポイントの寄与は、単一のポイントからそれを囲む空間領域に平滑化されます。個別に平滑化された寄与を集計すると、データの構造とその密度関数の全体像が得られます。以下では、このアプローチによって基礎となる密度関数の妥当な推定値が得られることを説明します。

前の図はカーネル密度推定のグラフ表現であり、ここではこれを正確に定義する。x ₁、x ₂、 ...、x _{n を}密度関数ƒで記述される共通分布から抽出されたd変量ランダムベクトルのサンプルとする。カーネル密度推定は次のように定義される。

${\displaystyle {\hat {f}}_{\mathbf {H} }(\mathbf {x} )={\frac {1}{n}}\sum _{i=1}^{n}K_{ \mathbf {H} }(\mathbf {x} -\mathbf {x} _{i})}$

どこ

• x = ( x ₁ , x ₂ , …, x _d ) ^T、 x _i = ( x _{i 1} , x _{i 2} , …, x _id ) ^T、i = 1, 2, …, nはdベクトルです
• Hは対称かつ正定値の帯域幅（または平滑化）d×d行列です。
• K は対称多変量密度であるカーネル関数です。
• ${\displaystyle K_{\mathbf {H} }(\mathbf {x} )=|\mathbf {H} |^{-1/2}K(\mathbf {H} ^{-1/2}\mathbf { x} )}$。

カーネル関数Kの選択はカーネル密度推定量の精度には重要ではないので、全体を通して標準の多変量正規カーネルを使用します。 、ここで H は共分散行列の役割を果たします。一方、帯域幅行列Hの選択は、誘導される平滑化の量と方向を制御するため、精度に影響を与える最も重要な要素です。^{[3] : 36–39}帯域幅行列も方向を誘導するという点は、1D カーネルでは方向が定義されていないため、多変量カーネル密度推定と単変量類似体との基本的な違いです。これにより、この帯域幅行列のパラメーター化が選択されます。3 つの主なパラメーター化クラス (複雑さの昇順) は、正のスカラーと単位行列の積のクラスS、主対角線上に正のエントリを持つ対角行列D 、対称正定値行列Fです。 Sクラスカーネルはすべての座標方向に同じ量のスムージングを適用し、Dカーネルは各座標で異なる量のスムージングを許可し、Fカーネルは任意の量と方向のスムージングを許可します。歴史的には、計算上の理由からSカーネルとDカーネルが最も普及していますが、研究によると、より一般的なFクラスカーネルを使用すると精度が大幅に向上することが示されています。^{[7] [8]}${\textstyle K_{\mathbf {H} }(\mathbf {x} )={(2\pi )^{-d/2}}\mathbf {|H|} ^{-1/2}e^{ -{\frac {1}{2}}\mathbf {x^{T}} \mathbf {H^{-1}} \mathbf {x} }}$

帯域幅行列を選択するための最も一般的に使用される最適基準は、MISEまたは平均積分二乗誤差である。

${\displaystyle \operatorname {MISE} (\mathbf {H} )=\operatorname {E} \!\left[\,\int ({\hat {f}}_{\mathbf {H} }(\mathbf {x} )-f(\mathbf {x} ))^{2}\,d\mathbf {x} \;\right].}$

これは一般に閉じた形式の表現を持たないので、その漸近近似（AMISE）を代理として使うのが普通である。

${\displaystyle \operatorname {AMISE} (\mathbf {H} )=n^{-1}|\mathbf {H} |^{-1/2}R(K)+{\tfrac {1}{4}}m_{2}(K)^{2}(\operatorname {vec} ^{T}\mathbf {H} )\mathbf {\Psi } _{4}(\operatorname {vec} \mathbf {H} )}$

どこ

• ${\displaystyle R(K)=\int K(\mathbf {x} )^{2}\,d\mathbf {x} }$、Kが正規核のとき、R ( K ) = (4 π ) ^{−d /2である。}
• ${\displaystyle \int \mathbf {x} \mathbf {x} ^{T}K(\mathbf {x} )\,d\mathbf {x} =m_{2}(K)\mathbf {I} _{d}}$、

ここで、I _{d は}d × dの 単位行列であり、正規カーネルではm _{2 = 1である。}

• D ² ƒ は、 ƒの2次偏微分のd × dヘッセ行列である。
• ${\displaystyle \mathbf {\Psi } _{4}=\int (\operatorname {vec} \,\operatorname {D} ^{2}f(\mathbf {x} ))(\operatorname {vec} ^{T}\operatorname {D} ^{2}f(\mathbf {x} ))\,d\mathbf {x} }$は、 ƒの4次偏微分を積分したd ² × d ²行列である。
• vecは行列の列を1つのベクトルに積み重ねるベクトル演算子です。例:${\displaystyle \operatorname {vec} {\begin{bmatrix}a&c\\b&d\end{bmatrix}}={\begin{bmatrix}a&b&c&d\end{bmatrix}}^{T}.}$

MISE ^{[3] : 97}に対するAMISE近似の品質は次のように与えられる。

${\displaystyle \operatorname {MISE} (\mathbf {H} )=\operatorname {AMISE} (\mathbf {H} )+o(n^{-1}|\mathbf {H} |^{-1/2}+\operatorname {tr} \,\mathbf {H} ^{2})}$

ここで、o は通常の小文字の o 表記を示します。経験的に、この記述は、サンプル サイズn → ∞ として AMISE が MISE の「良い」近似値であることを意味します。

任意の適切な帯域幅セレクタH は、要素ごとに大きな O 表記を適用してH = O ( n ^{−2/( d +4)} )であることが示されます。これを MISE 式に代入すると、最適な MISE はO ( n ^{−4/( d +4)} ) になります。^{[3] : 99–100}したがって、 n → ∞ では、MISE → 0、つまりカーネル密度推定値は平均二乗収束し、したがって確率的にも真の密度fに収束します。これらの収束モードは、カーネル法が適切な密度推定値につながるという動機のセクションでの記述を裏付けています。理想的な最適な帯域幅セレクタは

${\displaystyle \mathbf {H} _{\operatorname {AMISE} }=\operatorname {argmin} _{\mathbf {H} \in F}\,\operatorname {AMISE} (\mathbf {H} ).}$

この理想的なセレクターには未知の密度関数ƒが含まれているため、直接使用することはできません。AMISE のさまざまな推定値から、さまざまなデータ ベースの帯域幅セレクターが生まれます。ここでは、実際に最も広く適用できることがわかっている 2 つのクラスのセレクター、つまり平滑化クロス検証セレクターとプラグイン セレクターに焦点を当てます。

AMISEのプラグイン（PI）推定値は、Ψ _4をその推定値に置き換えることによって形成される。${\displaystyle {\hat {\mathbf {\Psi } }}_{4}}$

${\displaystyle \operatorname {PI} (\mathbf {H} )=n^{-1}|\mathbf {H} |^{-1/2}R(K)+{\tfrac {1}{4}}m_{2}(K)^{2}(\operatorname {vec} ^{T}\mathbf {H} ){\hat {\mathbf {\Psi } }}_{4}(\mathbf {G} )(\operatorname {vec} \,\mathbf {H} )}$

ここで である。したがって はプラグインセレクタである。^{[9] [10]これらの参考文献には、パイロット帯域幅行列}Gの最適推定アルゴリズムも含まれており、確率的にH _AMISEに収束することが示されている。 ${\displaystyle {\hat {\mathbf {\Psi } }}_{4}(\mathbf {G} )=n^{-2}\sum _{i=1}^{n}\sum _{j=1}^{n}[(\operatorname {vec} \,\operatorname {D} ^{2})(\operatorname {vec} ^{T}\operatorname {D} ^{2})]K_{\mathbf {G} }(\mathbf {X} _{i}-\mathbf {X} _{j})}$${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {PI} }=\operatorname {argmin} _{\mathbf {H} \in F}\,\operatorname {PI} (\mathbf {H} )}$${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {PI} }}$

平滑化クロスバリデーション（SCV）は、クロスバリデーション技術のより大きなクラスのサブセットです。SCV推定量は、プラグイン推定量と2番目の項で異なります。

${\displaystyle \operatorname {SCV} (\mathbf {H} )=n^{-1}|\mathbf {H} |^{-1/2}R(K)+n^{-2}\sum _{i=1}^{n}\sum _{j=1}^{n}(K_{2\mathbf {H} +2\mathbf {G} }-2K_{\mathbf {H} +2\mathbf {G} }+K_{2\mathbf {G} })(\mathbf {X} _{i}-\mathbf {X} _{j})}$

これがSCVセレクタです。^{[10] [11]} これらの参考文献には、パイロット帯域幅行列Gの最適推定アルゴリズムも含まれており、確率的にH _AMISEに収束することを確立しています。 ${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {SCV} }=\operatorname {argmin} _{\mathbf {H} \in F}\,\operatorname {SCV} (\mathbf {H} )}$${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {SCV} }}$

シルバーマンの経験則では、 を使用することが推奨されています。ここで、 は i 番目の変数の標準偏差、は次元数、 です。スコットの規則は です。 ${\displaystyle {\sqrt {\mathbf {H} _{ii}}}=\left({\frac {4}{d+2}}\right)^{\frac {1}{d+4}}n^{\frac {-1}{d+4}}\sigma _{i}}$${\displaystyle \sigma _{i}}$${\displaystyle d}$${\displaystyle \mathbf {H} _{ij}=0,i\neq j}$${\displaystyle {\sqrt {\mathbf {H} _{ii}}}=n^{\frac {-1}{d+4}}\sigma _{i}}$

最適帯域幅選択のセクションでは、MISEを導入しました。その構成は、密度推定量の期待値と分散に依存しています^{[3] ：97}

${\displaystyle \operatorname {E} {\hat {f}}(\mathbf {x} ;\mathbf {H} )=K_{\mathbf {H} }*f(\mathbf {x} )=f(\mathbf {x} )+{\frac {1}{2}}m_{2}(K)\operatorname {tr} (\mathbf {H} \operatorname {D} ^{2}f(\mathbf {x} ))+o(\operatorname {tr} \,\mathbf {H} )}$

ここで、*は2つの関数間の畳み込み演算子であり、

${\displaystyle \operatorname {Var} {\hat {f}}(\mathbf {x} ;\mathbf {H} )=n^{-1}|\mathbf {H} |^{-1/2}R(K)f(\mathbf {x} )+o(n^{-1}|\mathbf {H} |^{-1/2}).}$

これら2つの式が明確に定義されるためには、Hのすべての要素が0に近づくことと、nが無限大に近づくにつれてn ⁻¹ | H | ^−1/2が0に近づくことが必要である。これら2つの条件を仮定すると、期待値は真の密度fに近づく、つまりカーネル密度推定量は漸近的に不偏であり、分散は0に近づくことがわかる。標準的な平均二乗値分解を使用すると、

${\displaystyle \operatorname {MSE} \,{\hat {f}}(\mathbf {x} ;\mathbf {H} )=\operatorname {Var} {\hat {f}}(\mathbf {x} ;\mathbf {H} )+[\operatorname {E} {\hat {f}}(\mathbf {x} ;\mathbf {H} )-f(\mathbf {x} )]^{2}}$

MSE は 0 に近づくことがわかり、これはカーネル密度推定量が (平均二乗) 整合しており、したがって確率的に真の密度fに収束することを意味します。MSE の 0 への収束率は、前述の MISE 率O ( n ^{−4/(d+4) ) と必然的に同じであるため、密度推定量の}fへの収束率はO _p (n ^{−2/( d +4)} )です 。ここでO _{p は}確率の順序を表します。これにより、点ごとの収束が確立されます。関数の収束は、MISE の挙動を考慮し、十分な規則性があれば積分が収束率に影響を与えないことに注目することで同様に確立されます。

検討したデータに基づく帯域幅セレクタの場合、ターゲットはAMISE帯域幅マトリックスです。データに基づくセレクタは、相対速度O _p ( n ^{− α} )、α > 0 でAMISEセレクタに収束すると言います。

${\displaystyle \operatorname {vec} ({\hat {\mathbf {H} }}-\mathbf {H} _{\operatorname {AMISE} })=O(n^{-2\alpha })\operatorname {vec} \mathbf {H} _{\operatorname {AMISE} }.}$

プラグインセレクタと平滑化クロスバリデーションセレクタ（単一のパイロット帯域幅G_が与えられた場合）は両方とも相対速度Op（n ^{−2/（d +6）}）で収束することが確立されている^{[10] [12]}つまり、これらのデータベースセレクタは両方とも一貫した推定値です。

Rのksパッケージ^[13]は、プラグインと平滑化クロスバリデーションセレクター（その他）を実装しています。このデータセット（Rの基本配布に含まれています）には、米国イエローストーン国立公園のオールドフェイスフルガイザーの噴火継続時間（分）と次の噴火までの待ち時間（分）の2つの測定値を含む272のレコードが含まれています。

コード フラグメントは、プラグイン バンド幅マトリックスを使用してカーネル密度推定値を計算します。ここでも、色付きの等高線は、それぞれの確率質量を含む最小領域に対応します。赤 = 25%、オレンジ + 赤 = 50%、黄色 + オレンジ + 赤 = 75%。SCV セレクターを計算するには、を に置き換えます。これは、この例のプラグイン推定値とほぼ同じであるため、ここでは表示されません。 ${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {PI} }={\begin{bmatrix}0.052&0.510\\0.510&8.882\end{bmatrix}}.}$HpiHscv

ライブラリ( ks )
データ(忠実) 
H <- Hpi ( x =忠実) fhat <- kde ( x =忠実、H = H )プロット( fhat 、display = "filled.contour" 、drawpoints = TRUE 、cex = 0.5 、pch = 16 、col.pt = 1 )  
   
     

ランダム^に生成され_た500個_の点から、ガウス混合分布 (4π ) −1exp ⁽  −1⁄2 ( x12 + x22 ) ⁾ +( 4π ) ^−1exp  (−1⁄2 ( ( x1-3.5 ₎ 2 ⁺ x22 ) )の密度を推定することを検討する。2次元データにはMatlabルーチンを使用する。このルーチンは、2次ガウスカーネル用に特別に設計された自動帯域幅選択法である。[ ¹⁴ _]^図 は、自動的に選択された帯域幅を使用した結果の結合密度推定値を示している。

例の Matlab スクリプト

関数 kde2d.m を現在のディレクトリにダウンロードして保存した後、Matlab で次のコマンドを入力します。

  すべてクリア% 合成データを生成するdata =[ randn ( 500 , 2 ); randn ( 500 , 1 ) + 3.5 , randn ( 500 , 1 );]; % 現在のディレクトリに保存されているルーチンを呼び出す[帯域幅,密度, X , Y ]= kde2d ( data ); % データと密度推定値をプロットするcontour3 ( X , Y ,密度, 50 )、ホールドオンplot ( data (:, 1 ), data (:, 2 ), 'r.' , 'MarkerSize' , 5 )   
  
  
       
  
  
  
    
  

MISEは、密度推定値と真の密度関数fとの間の期待される積分L ₂距離である。これは、扱いやすさとほとんどのソフトウェアがMISEベースの帯域幅セレクタを実装していることから、最も広く使用されている。MISEが適切な尺度ではない場合をカバーする代替の最適性基準がある。^{[4] ：34–37、78}同等のL ₁尺度である平均積分絶対誤差は、

${\displaystyle \operatorname {MIAE} (\mathbf {H} )=\operatorname {E} \,\int |{\hat {f}}_{\mathbf {H} }(\mathbf {x} )-f(\mathbf {x} )|\,d\mathbf {x} .}$

その数学的解析はMISEのものよりかなり難しい。実際には、その利得は重要ではないようだ。^[15] L∞ノルムは平均一様絶対誤差である _。

${\displaystyle \operatorname {MUAE} (\mathbf {H} )=\operatorname {E} \,\operatorname {sup} _{\mathbf {x} }|{\hat {f}}_{\mathbf {H} }(\mathbf {x} )-f(\mathbf {x} )|.}$

^尤度誤差基準は、平均カルバック・ライブラー距離に基づくものを含む。

${\displaystyle \operatorname {MKL} (\mathbf {H} )=\int f(\mathbf {x} )\,\operatorname {log} [f(\mathbf {x} )]\,d\mathbf {x} -\operatorname {E} \int f(\mathbf {x} )\,\operatorname {log} [{\hat {f}}(\mathbf {x} ;\mathbf {H} )]\,d\mathbf {x} }$

平均ヘリンガー距離

${\displaystyle \operatorname {MH} (\mathbf {H} )=\operatorname {E} \int ({\hat {f}}_{\mathbf {H} }(\mathbf {x} )^{1/2}-f(\mathbf {x} )^{1/2})^{2}\,d\mathbf {x} .}$

KLはクロスバリデーション法を用いて推定することができるが、KLクロスバリデーションセレクターは、たとえ有界密度関数に対して一貫性を保っていたとしても、最適ではない可能性がある。 ^[17] MHセレクターは文献で簡単に検討されている。^[18]

これらの最適性基準はすべて距離に基づいた尺度であり、必ずしもより直感的な近さの概念に対応するわけではないため、この懸念に応じてより視覚的な基準が開発されてきました。^[19]

最近の研究では、カーネルとその帯域幅は、分布の形について何の仮定もせずに、入力データ自体から最適かつ客観的に選択できることが示されています。^[20] 結果として得られるカーネル密度推定値は、サンプルが追加されるにつれて、真の確率分布に急速に収束します。これは、パラメトリック推定値に期待される速度に近いものです。^{[20] [21] [22]} このカーネル推定値は、単変量サンプルでも多変量サンプルでも同様に機能します。最適なカーネルは、フーリエ空間で、データのフーリエ変換、経験的特性関数の観点から、最適な減衰関数（カーネルのフーリエ変換）として定義されます（カーネル密度推定を参照）。 ${\displaystyle n^{-1}}$${\displaystyle {\hat {\psi _{h}}}({\vec {t}})}$${\displaystyle {\hat {K}}({\vec {x}})}$${\displaystyle {\hat {\varphi }}({\vec {t}})}$

${\displaystyle {\hat {\psi _{h}}}({\vec {t}})\equiv {\frac {N}{2(N-1)}}\left[1+{\sqrt {1-{\frac {4(N-1)}{N^{2}|{\hat {\varphi }}({\vec {t}})|^{2}}}}}I_{\vec {A}}({\vec {t}})\right]}$ ^[22]

${\displaystyle {\hat {f}}(x)={\frac {1}{(2\pi )^{d}}}\int {\hat {\varphi }}({\vec {t}})\psi _{h}({\vec {t}})e^{-i{\vec {t}}\cdot {\vec {x}}}d{\vec {t}}}$

ここで、Nはデータ ポイントの数、dは次元 (変数) の数、は「受け入れられる頻度」の場合は 1、それ以外の場合は 0 に等しいフィルターです。このフィルター関数を定義する方法はさまざまですが、単変量または多変量サンプルに機能する単純な方法は「最低連続ハイパーボリューム フィルター」と呼ばれ、受け入れられる頻度が原点を囲む連続した周波数のサブセットのみになるように選択されます(このフィルター関数と他のフィルター関数については ^{[22]を参照してください)。}${\displaystyle I_{\vec {A}}({\vec {t}})}$${\displaystyle I_{\vec {A}}({\vec {t}})}$${\displaystyle |{\hat {\varphi }}({\vec {t}})|^{2}\geq 4(N-1)N^{-2}}$

経験的特性関数（ECF）の直接計算は、本質的にデータサンプルの直接フーリエ変換を伴うため、遅いことに注意してください。しかし、ECFは非均一高速フーリエ変換（nuFFT）法を使用して正確に近似できることがわかっています。 ^{[21] [22]}これにより、計算速度が数桁向上します（問題の次元によって異なります）。この客観的KDE法とnuFFTベースのECF近似の組み合わせは、文献ではfastKDEと呼ばれています。 ^[22]

• カーネル密度推定 – 単変量カーネル密度推定。
• 可変カーネル密度推定 – 可変帯域幅のカーネルを使用した多変量密度の推定

• 「Multivariate Kernel Smoothing and Its Applications」は、密度推定を含むカーネルスムージングの多くのトピックを網羅した包括的な書籍です。R の ks パッケージ コード スニペットが含まれています。
• kde2d.m 二変量カーネル密度推定のためのMatlab関数。
• libagf多変量、可変帯域幅カーネル密度推定用のC++ライブラリ。
• akde.m多変量、可変帯域幅カーネル密度推定用のMatlab m ファイル。
• PyQt-Fit パッケージの helit および pyqt_fit.kde モジュールは、多変量カーネル密度推定用のPythonライブラリです。