カーネル密度推定

統計学において、カーネル密度推定( KDE ) は、カーネル平滑化を確率密度推定に適用したもの、つまり、重みとしてカーネルに基づいてランダム変数の確率密度関数を推定するノンパラメトリックな方法です。KDE は、有限のデータサンプルに基づいて母集団についての推論が行われる基本的なデータ平滑化問題に答えます。信号処理や計量経済学などの分野では、エマニュエル・パルゼンとマレー・ローゼンブラットにちなんで、パルゼン・ローゼンブラットウィンドウ法とも呼ばれ、通常、この 2 人が現在の形式を独立に考案したとされています。^{[1] [2]カーネル密度推定の有名な応用例の 1 つは、}ナイーブベイズ分類器を使用するときにデータのクラス条件付き周辺密度を推定することであり、これにより予測精度が向上します。^[3]

( x ₁ , x ₂ , ..., x _n )を、任意の点xにおける未知の密度ƒを持つ一変量分布から抽出された、独立かつ同一に分布するサンプルとします。この関数ƒの形状を推定することに関心があります。そのカーネル密度推定量は

${\displaystyle {\widehat {f}}_{h}(x)={\frac {1}{n}}\sum _{i=1}^{n}K_{h}(x-x_{i})={\frac {1}{nh}}\sum _{i=1}^{n}K{\Big (}{\frac {x-x_{i}}{h}}{\Big )},}$

ここで、Kはカーネル（非負関数）であり、 h > 0は帯域幅または単に幅と呼ばれる平滑化パラメータです。[3] 下付き文字hを持つカーネルはスケールカーネルと呼ばれ、K h ( x ) = K ( ) と定義されます。直感^的には、データが_許す限り小さいhを選択したいものですが、推定量のバイアスと分散の間には常にトレードオフが存在します。帯域幅の選択については、以下で詳しく説明します。 ${\displaystyle {\tfrac {1}{h}}}$ ${\displaystyle {\tfrac {x}{h}}}$

カーネル関数には、均一、三角、バイウェイト、トリウェイト、エパネチニコフ（放物線）、正規など、さまざまな関数が一般的に使用されています。エパネチニコフカーネルは平均二乗誤差の意味で最適ですが、[ ^4]前述のカーネルでは効率の低下は小さいです。^[5]数学的な特性が便利なため、正規カーネルがよく使用されます。これは、 K ( x ) = ϕ ( x )を意味します。ここで、ϕ は標準正規密度関数です。カーネル密度推定量は次のようになります。

${\displaystyle {\widehat {f}}_{h}(x)={\frac {1}{nh\sigma }}{\frac {1}{\sqrt {2\pi }}}\sum _{i=1}^{n}\exp {\Big (}{\frac {-(x-x_{i})^{2}}{2h^{2}\sigma ^{2}}}{\Big )},}$

ここで、サンプルの標準偏差です。 ${\displaystyle \sigma}$${\displaystyle {\vec {x}}}$

カーネル密度推定の構築は、密度推定以外の分野でも解釈が見出されます。^[6]たとえば、熱力学では、これは熱カーネル（熱方程式の基本解）が各データポイントの位置x _iに配置されたときに発生する熱量に相当します。同様の方法は、多様体学習（拡散マップなど）のための点群上に離散ラプラス演算子を構築するのにも使用されます。

カーネル密度推定はヒストグラムと密接に関連していますが、適切なカーネルを使用することで、滑らかさや連続性などの特性を付与することができます。以下の図は、これら 6 つのデータ ポイントに基づいてこの関係を示しています。

ヒストグラムの場合、まず、水平軸がデータの範囲をカバーするサブ間隔またはビンに分割されます。この場合、幅 2 のビンが 6 つあります。データ ポイントがこの間隔内に入ると、高さ 1/12 のボックスがそこに配置されます。複数のデータ ポイントが同じビン内に入る場合、ボックスは互いに積み重ねられます。

カーネル密度推定では、標準偏差1.5の正規カーネル（赤い破線で示される）が各データポイントx _i上に配置されます。カーネルを合計してカーネル密度推定値（青い実線）が作成されます。カーネル密度推定値の滑らかさ（ヒストグラムの離散性と比較）は、連続ランダム変数の場合、カーネル密度推定値が真の基礎密度にどれだけ速く収束するかを示しています。^[7]

カーネルの帯域幅は、結果の推定値に強い影響を与える自由パラメータです。その影響を示すために、標準正規分布からシミュレートされたランダム サンプルを取得します(横軸のラグ プロットの青いスパイクにプロットされています)。灰色の曲線は、真の密度 (平均 0、分散 1 の正規密度) です。これに比べて、赤い曲線は、帯域幅h = 0.05 という小さすぎる帯域幅を使用したために生じた偽のデータ アーティファクトが多すぎるため、平滑化が不十分です。緑の曲線は、帯域幅h = 2を使用すると基礎構造の多くが不明瞭になるため、平滑化が過剰です。帯域幅h = 0.337 の黒い曲線は、密度推定値が真の密度に近いため、最適に平滑化されていると考えられます。極端な状況は、限界(平滑化なし) で発生します。この限界では、推定値は、分析されたサンプルの座標を中心とするn個のデルタ関数の合計になります。他の極端な限界では、推定値はサンプルの平均を中心とした、使用されたカーネルの形状を保持します（完全に滑らか）。 ${\displaystyle h\to 0}$ ${\displaystyle h\to \infty}$

このパラメータを選択するために使用される最も一般的な最適性基準は、期待されるL ₂ リスク関数（平均積分二乗誤差とも呼ばれる）です。

${\displaystyle \operatorname {MISE} (h)=\operatorname {E} \!\left[\,\int ({\hat {f}}_{h}(x)-f(x))^{2}\,dx\right]}$

ƒ とKに関する弱い仮定の下では（ƒ は一般には知られていない実密度関数である）、^{[1] [2]}

${\displaystyle \operatorname {MISE} (h)=\operatorname {AMISE} (h)+{\mathcal {o}}((nh)^{-1}+h^{4})}$

ここでoは小文字のo表記、nはサンプル数（上記と同じ）である。AMISEは漸近的MISE、すなわち2つの主要な項である。

${\displaystyle \operatorname {AMISE} (h)={\frac {R(K)}{nh}}+{\frac {1}{4}}m_{2}(K)^{2}h^{4}R(f'')}$

ここで関数gに対して、 は関数の2次導関数であり、は核である。このAMISEの最小値は、この微分方程式の解である。 ${\displaystyle R(g)=\int g(x)^{2}\,dx}$${\displaystyle m_{2}(K)=\int x^{2}K(x)\,dx}$${\displaystyle f''}$${\displaystyle f}$${\displaystyle K}$

${\displaystyle {\frac {\partial }{\partial h}}\operatorname {AMISE} (h)=-{\frac {R(K)}{nh^{2}}}+m_{2}(K)^{2}h^{3}R(f'')=0}$

または

${\displaystyle h_{\operatorname {AMISE} }={\frac {R(K)^{1/5}}{m_{2}(K)^{2/5}R(f'')^{1/5}}}n^{-1/5}=Cn^{-1/5}}$

AMISE式もh _AMISE式も、未知の密度関数またはその2次導関数を含むため、直接使用することはできません。この困難を克服するために、帯域幅を選択するためのさまざまな自動データベース手法が開発されました。それらの有効性を比較するレビュー研究がいくつか行われており、^{[8] [9] [10] [11] [12] [13] [14]}、プラグインセレクター^{[6] [15] [16]}とクロスバリデーションセレクター^{[17] [18] [19]}が幅広いデータセットで最も有用であるという一般的なコンセンサスが得られています。 ${\displaystyle f}$${\displaystyle f''}$

_AMISEのhと同じ漸近順序n ^−1/5を持つ任意の帯域幅hをAMISEに代入すると、AMISE( h ) = O ( n ^−4/5 )が得られます。ここで、Oは大きなO表記です。弱い仮定の下では、カーネル推定量よりも速い速度で収束するノンパラメトリック推定量は存在し得ないことが示されます。^[20] n ^{−4/5の収束率は、パラメトリック法の典型的な}n ⁻¹収束率よりも遅いことに注意してください。

帯域幅が固定されず、推定値（バルーン推定器）またはサンプル（ポイントワイズ推定器）の位置に応じて変化する場合は、適応型または可変帯域幅カーネル密度推定と呼ばれる特に強力な方法が生成されます。

重い裾を持つ分布のカーネル密度推定のためのバンド幅の選択は比較的困難である。^[21]

ガウス基底関数が単変量データの近似に使用され、推定される基礎密度がガウスである場合、h（つまり、平均積分二乗誤差を最小化する帯域幅）の最適な選択は次のようになる。^[22]

${\displaystyle h=\left({\frac {4{\hat {\sigma }}^{5}}{3n}}\right)^{\frac {1}{5}}\approx 1.06\,{\hat {\sigma }}\,n^{-1/5},}$

値がより堅牢であるとは、ロングテール分布や歪んだ分布、または二峰性混合分布への適合性が向上する場合です。これは、標準偏差を以下の パラメータに置き換えることによって経験的に行われることが多いです。${\displaystyle h}$ ${\displaystyle {\hat {\sigma }}}$${\displaystyle A}$

${\displaystyle A=\min \left({\hat {\sigma }},{\frac {IQR}{1.34}}\right)}$ここで、IQR は四分位範囲です。

モデルを改善するもう 1 つの変更は、係数を 1.06 から 0.9 に減らすことです。最終的な式は次のようになります。

${\displaystyle h=0.9\,\min \left({\hat {\sigma }},{\frac {IQR}{1.34}}\right)\,n^{-{\frac {1}{5}}}}$

サンプルサイズは どこですか。${\displaystyle n}$

この近似は正規分布近似、ガウス近似、またはシルバーマンの経験則と呼ばれます。^[22]この経験則は計算が簡単ですが、密度が正規分布に近くない場合に大幅に不正確な推定値をもたらす可能性があるため、注意して使用する必要があります。たとえば、二峰性ガウス混合モデルを推定する場合

${\displaystyle {\frac {1}{2{\sqrt {2\pi }}}}e^{-{\frac {1}{2}}(x-10)^{2}}+{\frac {1}{2{\sqrt {2\pi }}}}e^{-{\frac {1}{2}}(x+10)^{2}}}$

右の図は、200点のサンプルから得られた真の密度と2つのカーネル密度推定値を示しています。1つは経験則帯域幅を使用し、もう1つは方程式を解く帯域幅を使用しています。^{[6] [16]}経験則帯域幅に基づく推定値は、大幅に平滑化されています。

サンプル（x ₁、x ₂、...、x _n ）が与えられれば、特性関数 φ（t）= E[ e ^itX ]を次のように 推定するのが自然である。

${\displaystyle {\widehat {\varphi }}(t)={\frac {1}{n}}\sum _{j=1}^{n}e^{itx_{j}}}$

特性関数がわかれば、フーリエ変換式によって対応する確率密度関数を求めることができます。この逆変換式を適用する際の難しさの 1 つは、大きなtに対して推定値が信頼できないため、積分が発散してしまうことです。この問題を回避するには、推定値に減衰関数ψ _h ( t ) = ψ ( ht )を乗算します。この関数は原点で 1 になり、無限大で 0 になります。「帯域幅パラメーター」h は、関数を減衰させる速度を制御します。特にhが小さい場合、ψ _h ( t ) はtの広い範囲でほぼ 1 になります。つまり、tの最も重要な領域では実質的に変化しません。 ${\displaystyle \scriptstyle {\widehat {\varphi }}(t)}$${\displaystyle \scriptstyle {\widehat {\varphi }}(t)}$${\displaystyle \scriptstyle {\widehat {\varphi }}(t)}$${\displaystyle \scriptstyle {\widehat {\varphi }}(t)}$

関数ψの最も一般的な選択は、一様関数ψ ( t ) = 1 {−1 ≤ t ≤ 1 } （これは、逆関数の積分区間を[−1/ h , 1/ h ]に切り捨てることを意味します）またはガウス関数 ψ ( t ) = e ^{− π t 2 の}いずれかです。関数ψが選択されると、逆関数の式が適用され、密度推定量は次のように表されます 。

${\displaystyle {\begin{aligned}{\widehat {f}}(x)&={\frac {1}{2\pi }}\int _{-\infty }^{+\infty }{\widehat {\varphi }}(t)\psi _{h}(t)e^{-itx}\,dt={\frac {1}{2\pi }}\int _{-\infty }^{+\infty }{\frac {1}{n}}\sum _{j=1}^{n}e^{it(x_{j}-x)}\psi (ht)\,dt\\[5pt]&={\frac {1}{nh}}\sum _{j=1}^{n}{\frac {1}{2\pi }}\int _{-\infty }^{+\infty }e^{-i(ht){\frac {x-x_{j}}{h}}}\psi (ht)\,d(ht)={\frac {1}{nh}}\sum _{j=1}^{n}K{\Big (}{\frac {x-x_{j}}{h}}{\Big )},\end{aligned}}}$

ここで、Kは減衰関数ψのフーリエ変換です。したがって、カーネル密度推定量は特性関数密度推定量と一致します。

(グローバル) モードの定義をローカルな意味に拡張し、ローカル モードを定義することができます。

${\displaystyle M=\{x:g(x)=0,\lambda _{1}(x)<0\}}$

つまり、は密度関数が局所的に最大になる点の集合です。 の自然な推定量はKDEのプラグインです^{[23] [24]。}ここで、およびはKDE版のおよびです。 軽度の仮定の下では、は矛盾のない推定量です。 平均シフトアルゴリズム^{[25] [26] [27]}を使用して推定量を数値的に計算できることに留意してください。 ${\displaystyle M}$${\displaystyle M}$${\displaystyle g(x)}$${\displaystyle \lambda _{1}(x)}$${\displaystyle g(x)}$${\displaystyle \lambda _{1}(x)}$${\displaystyle M_{c}}$${\displaystyle M}$${\displaystyle M_{c}}$

カーネル密度推定器のソフトウェア実装の非網羅的なリストには以下が含まれます。

• Analyticaリリース 4.4では、PDF 結果のスムージングオプションは KDE を使用し、式からは組み込み関数を介して利用できますPdf。
• C / C++では、FIGTree は通常のカーネルを使用してカーネル密度推定値を計算するために使用できるライブラリです。MATLAB インターフェイスが利用可能です。
• C++では、 libagf は可変カーネル密度推定用のライブラリです。
• C++では、mlpack はさまざまなカーネルを使用して KDE を計算できるライブラリです。これにより、計算を高速化するために誤差許容値を設定できます。PythonおよびRインターフェイスが利用可能です。
• C#とF#では、Math.NET Numericsはカーネル密度推定を含む数値計算のためのオープンソースライブラリです。
• CrimeStatでは、カーネル密度推定は、正規、一様、4 次、負の指数、および三角の 5 つの異なるカーネル関数を使用して実装されています。シングルカーネル密度推定ルーチンとデュアルカーネル密度推定ルーチンの両方が利用可能です。カーネル密度推定は、Head Bang ルーチンの補間、2 次元の Journey-to-crime 密度関数の推定、および 3 次元の Bayesian Journey-to-crime 推定の推定にも使用されます。
• ELKIでは、カーネル密度関数はパッケージ内にあります。de.lmu.ifi.dbs.elki.math.statistics.kernelfunctions
• ESRI製品では、カーネル密度マッピングは Spatial Analyst ツールボックスの外部で管理され、Quartic(biweight) カーネルが使用されます。
• Excelでは、王立化学協会が分析方法委員会の技術概要 4 に基づいてカーネル密度推定を実行するアドインを作成しました。
• gnuplotでは、カーネル密度推定はsmooth kdensityオプションによって実装されており、データファイルには各ポイントの重みと帯域幅を含めることができ、または帯域幅は「シルバーマンの経験則」に従って自動的に設定することもできます^[28]（上記参照）。
• Haskellでは、カーネル密度は統計パッケージに実装されています。
• IGOR Proでは、カーネル密度推定は操作によって実装されます(Igor Pro 7.00 で追加)。帯域幅はユーザーが指定するか、Silverman、Scott、または Bowmann とAzzaliniStatsKDEによって推定できます。カーネル タイプは、Epanechnikov、Bi-weight、Tri-weight、Triangular、Gaussian、および Rectangular です。
• Javaでは、Weka機械学習パッケージは、weka.estimators.KernelEstimator などを提供します。
• JavaScriptでは、視覚化パッケージD3.js がscience.stats パッケージ内に KDE パッケージを提供しています。
• JMPのグラフビルダー プラットフォームでは、カーネル密度推定を利用して、二変量密度の等高線図と高密度領域 (HDR) と、一変量密度のバイオリン プロットと HDR を提供します。スライダーを使用すると、帯域幅を変更できます。二変量と一変量のカーネル密度推定は、それぞれ Y と X の相関と分布のプラットフォームでも提供されます。
• Juliaでは、カーネル密度推定は KernelDensity.jl パッケージに実装されています。
• KNIMEでは、1D カーネル密度プロットなどのVernalisコミュニティ貢献のノードを使用して、1D および 2D カーネル密度分布を生成およびプロットできます。基礎となる実装はJavaで記述されています。
• MATLABでは、カーネル密度推定はksdensity関数（統計ツールボックス）を通じて実装されています。MATLABの2018aリリースでは、帯域幅とカーネルスムーザーの両方を指定でき、カーネル密度の範囲を指定するなどの他のオプションも指定できます。^[29]また、自動帯域幅選択法を実装した無料のMATLABソフトウェアパッケージ^[6]は、MATLAB Central File Exchangeから入手できます。
  • 1次元データ
  • 2次元データ
  • n次元データカーネル回帰、カーネル密度推定、ハザード関数のカーネル推定などの実装を備えた無料のMATLABツールボックスがこれらのページで入手できます（このツールボックスは書籍^[30]
    の一部です）。
• Mathematicaでは、数値カーネル密度推定は関数SmoothKernelDistribution^[31]によって実装され、記号推定は関数KernelMixtureDistribution^[32]を使用して実装されており、どちらもデータ駆動型の帯域幅を提供します。
• Minitabでは、王立化学協会が分析方法委員会の技術概要4に基づいてカーネル密度推定を実行するマクロを作成しました。^[33]
• NAGライブラリでは、カーネル密度推定はルーチン（ライブラリのg10baFortran ^[34]とC ^{[​​35] の両方のバージョンで利用可能）を介して実装されています。}
• Nuklei では、C++カーネル密度法は特殊ユークリッド群のデータに焦点を当てます。${\displaystyle SE(3)}$
• Octaveでは、カーネル密度推定はkernel_densityオプション (計量経済学パッケージ) によって実装されます。
• Originでは、ユーザー インターフェイスから 2D カーネル密度プロットを作成でき、LabTalk、 Python、またはCコードから 1D 用の Ksdensity と 2D 用の Ks2density という 2 つの関数を使用できます。
• Perlでは、Statistics-KernelEstimationモジュールに実装があります。
• PHPでは、実装はMathPHPライブラリにあります。
• Pythonには、多くの実装が存在します。PyQt-Fit パッケージの pyqt_fit.kde モジュール、SciPy ( scipy.stats.gaussian_kde)、Statsmodels (KDEUnivariateおよびKDEMultivariate)、scikit-learn ( KernelDensity) などです (比較を参照^[36] )。KDEpy は重み付きデータをサポートし、その FFT 実装は他の実装よりも桁違いに高速です。一般的に使用されている pandas ライブラリ [1] は、plot メソッド ( df.plot(kind='kde')[2] ) を通じて kde プロットをサポートしています。重み付き相関 MCMC サンプル用の getdist パッケージは、1D および 2D 分布の最適化された帯域幅、境界補正、および高次法をサポートしています。カーネル密度推定用の新しく使用されたパッケージの 1 つは seaborn ( import seaborn as sns、sns.kdeplot()) です。^[37] KDE の GPU 実装も存在します。^[38]
• Rでは、density基本分布でを通じて実装され、bw.nrd0関数は stats パッケージで使用されます。この関数は、Silverman の書籍にある最適化された式を使用します。bkdeKernSmooth ライブラリ、 ParetoDensityEstimationDataVisualizations ライブラリ (パレート分布密度推定用)、kdeks ライブラリ、dkdenevmixdbckdenライブラリ (後者は境界サポートの境界補正カーネル密度推定用)、npudensnp ライブラリ (数値およびカテゴリデータ)、sm.densitysm ライブラリで使用されます。 パッケージやライブラリのインストールを必要としない関数の実装についてはkde.R、kde.R を参照してください。都市分析専用の btb ライブラリは、を通じてカーネル密度推定を実装しますkernel_smoothing。
• SASでは、proc kde単変量および二変量カーネル密度を推定するために使用できます。
• Apache Sparkでは、KernelDensity()クラス^[39]
• Stataでは、例えば^[40]kdensityを通じて実装されています。あるいは、無料のStataモジュールKDENS ^[41]を使用して、ユーザーが1Dまたは2Dの密度関数を推定できるようにすることもできます。histogram x, kdensity
• Swiftでは、SwiftStats.KernelDensityEstimationオープンソースの統計ライブラリ SwiftStats を通じて実装されています。

ウィキメディア・コモンズには、カーネル密度推定に関連するメディアがあります。

• カーネル（統計）
• カーネルスムージング
• カーネル回帰
• 密度推定（他の例の提示を含む）
• 平均シフト
• スケール空間: 3 つ組 {( x、h 、 xで評価された帯域幅hを持つ KDE : すべてx、h > 0}は、データのスケール空間表現を形成します。
• 多変量カーネル密度推定
• 可変カーネル密度推定
• ヘッド/テールブレイク

• Härdle、Müller、Sperlich、Werwatz、「ノンパラメトリックおよびセミパラメトリック手法」、Springer-Verlag Berlin Heidelberg 2004、pp. 39–83

• カーネル密度推定の紹介 ヒストグラムの改良としてカーネル密度推定器を説明する短いチュートリアルです。
• カーネル帯域幅の最適化 最適化されたカーネル密度の推定値を生成する無料のオンライン ツール。
• 無料のオンライン ソフトウェア (計算機) は、ガウス、エパネチニコフ、長方形、三角形、バイウェイト、コサイン、オプトコサインの各カーネルに従って、データ シリーズのカーネル密度推定を計算します。
• カーネル密度推定アプレット カーネル密度推定のオンライン インタラクティブな例。.NET 3.0 以降が必要です。