手元に身長のデータが100人分あるとします。「このデータはどんな分布から来たのか」を知りたい。いちばん手軽なのはヒストグラムを描くことです。
ところが描いてみると困ったことが起きます。棒の幅を変えると形がガラッと変わる。棒の境目をどこに置くかでも見え方が違う。そのうえ、本来なめらかなはずの分布が、階段みたいなギザギザになってしまう。「データの裏にある、なめらかな確率密度」をそのまま推定する方法はないものでしょうか。
そこで登場するのがカーネル密度推定（Kernel Density Estimation, KDE）です。考え方はとても素直で、「各データ点の上に小さな山をひとつずつ置き、それを全部足し合わせて1本のなめらかな曲線にする」だけ。まずはこのイメージを絵で掴んでください。
黒い縦棒がデータ点です。その1つひとつの上に、青い山（カーネル）を置きます。点が近くに集まっている場所では山が重なって高くなり、まばらな場所では低いまま。これを全部足したのがオレンジの曲線、つまり密度の推定結果です。ヒストグラムのように「箱に放り込んで数える」のではなく、「各点のまわりに影響をなめらかに広げて重ねる」。これだけで、棒の境目の恣意性も階段状の段差も消えてしまいます。母集団の分布の形をあらかじめ仮定しない、ノンパラメトリックな手法の代表格です。
カーネル密度推定がわかると、こんな場面で役に立ちます。
- データの可視化: ヒストグラムより安定して分布の形をつかめる。violin plot や seaborn の
kdeplotの中身はまさに KDE です。 - 異常検知: 正常データの密度を KDE で推定し、密度が低い場所に落ちる点を「異常」とみなす（密度ベース異常検知の基礎）。
- ノンパラメトリックなベイズ・生成: クラスごとの密度を KDE で推定するナイーブベイズの連続版や、データ分布からのリサンプリングに使えます。
本記事で扱う内容は次のとおりです。
- ヒストグラムの欠点を出発点に、KDE の定義を直感から導く
- カーネル関数とバンド幅の役割、それぞれが満たすべき条件を整理する
- KDE 推定量のバイアスと分散を1行ずつ導出し、AMISE（漸近平均積分二乗誤差）にまとめる
- AMISE を最小にする最適バンド幅と、実務で使う Silverman の目安公式を導く
- Python で KDE をスクラッチ実装し、図とともにバンド幅の効果・2次元 KDE・scipy との一致を確かめる
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 確率変数と確率密度関数 — 推定したい対象である確率密度の定義
- 正規分布とは？定義やグラフの見方を徹底解説 — いちばんよく使うガウスカーネルの土台
- バイアス・バリアンスのトレードオフ（統計） — バンド幅選びの核心はこのトレードオフ
- 混合分布とは？混合正規分布の理論とEM推定をわかりやすく解説 — KDE は「n個のカーネルを等しい重みで混ぜた分布」とも読めます
テイラー展開、期待値・分散の計算、簡単な置換積分ができれば、導出は最後まで追えます。
ヒストグラムからカーネル密度推定へ
まず、推定したい対象をはっきりさせましょう。観測データ $x_1, x_2, \dots, x_n$ が、ある未知の確率密度関数 $f(x)$ から独立に得られたとします。ゴールは、このデータだけから $f(x)$ を復元することです。
いちばん古典的な推定はヒストグラムです。区間を幅 $h$ のビンに区切り、各ビンに入った点の数を数え、全体の積分が1になるように正規化する。確かに密度の推定にはなりますが、2つの不満が残ります。ひとつは、ビンの境界の位置を少しずらすだけで、同じ点が隣のビンに移って形が変わってしまうこと。もうひとつは、推定値がビンごとに一定なので、密度が階段状に不連続になってしまうことです。
左と中央は、まったく同じデータに対してビンの境界を $0.5$ だけずらしたヒストグラムです。山の高さも谷の位置も変わってしまい、「どちらが本当の分布か」を決められません。右の KDE は境界という概念を持たないので、こうしたブレがなく、なめらかな1本の曲線に落ち着きます。この安定性こそ KDE を使う最大の動機です。
不連続の原因は「点をビンという固定の箱に放り込む」ところにあります。発想を変えましょう。箱を固定するのではなく、各データ点を中心に小さな山をひとつずつ置き、それを全部足し合わせたらどうなるか。点が密集する場所では山が何重にも重なって高くなり、まばらな場所では低いまま。これなら境界の恣意性はないし、山をなめらかな関数にしておけば結果もなめらかになります。冒頭の絵で見たとおりです。
この「各点に置く山」をカーネル関数 $K(\cdot)$ と呼び、山の幅をバンド幅 $h$ で決めます。すると密度推定量は次のように書けます。
$$ \begin{equation} \hat{f}_h(x) = \frac{1}{nh} \sum_{i=1}^{n} K\!\left(\frac{x – x_i}{h}\right) \end{equation} $$
各項 $K((x-x_i)/h)$ が「データ点 $x_i$ を中心に幅 $h$ で広がる山」、それを $n$ 点ぶん足して $nh$ で割る。$1/(nh)$ は積分を1に保つための正規化です。次は、この式の部品であるカーネルとバンド幅が何を満たすべきかを整理します。
カーネル関数とバンド幅
カーネル関数の条件
$\hat{f}_h(x)$ がまっとうな確率密度（非負で積分1）になるために、カーネル $K$ にも条件を課します。標準的には、$K$ 自身が平均0・分散有限の対称な確率密度であることを求めます。式で書くと次の3つです。
$$ \int K(u)\,du = 1, \qquad \int u\,K(u)\,du = 0, \qquad \int u^2 K(u)\,du = \mu_2(K) < \infty $$
1番目は $\hat{f}_h$ の積分が1になることを保証し、2番目（対称性）は推定のバイアスを小さく保ち、3番目の $\mu_2(K)$ は後でバイアスの大きさに効いてきます。よく使うカーネルは次のとおりです。
- ガウスカーネル: $K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}$。なめらかで裾が無限に伸び、いちばん標準的。
- Epanechnikov カーネル: $K(u) = \frac{3}{4}(1-u^2)\,\mathbb{1}[|u|\le 1]$。後で出る AMISE を最小にする「最適なカーネル」として知られます。
- 一様（箱型）カーネル: $K(u) = \frac{1}{2}\mathbb{1}[|u|\le 1]$。これを使うとヒストグラムに近くなります（ただし境界の恣意性はない）。
3つのカーネルの形を並べました。どれも「中央が高く、左右対称で、面積が1」という共通の性質を持ち、違うのは裾の広がり方と滑らかさだけです。ガウスは裾が無限に伸びてなめらか、Epanechnikov と一様は $|u|>1$ で完全に0になります。面白いことに、推定精度はカーネルの種類にはほとんど左右されません。本当に効くのは次のバンド幅です。
バンド幅の役割
バンド幅 $h$ は山の幅、つまり「各点の影響をどれだけ広げるか」を決めます。これが KDE で最も重要なパラメータです。
- $h$ が小さすぎると、各点が針のように鋭い山になり、推定はデータの偶然のばらつきまで拾ってギザギザになります（過小平滑化、分散が大きい）。
- $h$ が大きすぎると、山が広がりすぎて互いに溶け合い、本来あるはずの山や谷がのっぺり潰れます（過剰平滑化、バイアスが大きい）。
つまりバンド幅選びは、典型的なバイアスと分散のトレードオフです。では「ちょうどよい $h$」を理論的に決められるのか。それには、推定量 $\hat{f}_h(x)$ が真の $f(x)$ からどれだけずれるかを定量化する必要があります。次の節でバイアスと分散を導出します。
推定量のバイアスと分散を導出する
KDE の良し悪しを測るため、各点 $x$ での推定量 $\hat{f}_h(x)$ の期待値と分散を、真の密度 $f$ で表しましょう。$x_1,\dots,x_n$ は $f$ からの独立標本なので、$\hat{f}_h(x)$ は確率変数であり、そのばらつきを評価します。
バイアスの導出
まず期待値を計算します。$\hat{f}_h(x)$ は $n$ 個の同分布な項の平均なので、期待値は1項の期待値に等しくなります。
$$ \mathbb{E}[\hat{f}_h(x)] = \frac{1}{h}\,\mathbb{E}\!\left[K\!\left(\frac{x – X}{h}\right)\right] = \frac{1}{h}\int K\!\left(\frac{x – t}{h}\right) f(t)\,dt $$
ここで $X$ は $f$ に従う1つの観測を表し、期待値を密度 $f(t)$ による積分で書き下しました。この積分を扱いやすくするため、置換 $u = (t – x)/h$、すなわち $t = x + hu$、$dt = h\,du$ を行います。すると $1/h$ と $dt$ の $h$ が打ち消し合い、
$$ \mathbb{E}[\hat{f}_h(x)] = \int K(-u)\, f(x + hu)\,du = \int K(u)\, f(x + hu)\,du $$
となります（最後にカーネルの対称性 $K(-u)=K(u)$ を使いました）。次に、$h$ が小さいことを使って $f(x+hu)$ を $x$ のまわりでテイラー展開します。
$$ f(x + hu) = f(x) + hu\,f'(x) + \frac{1}{2}h^2 u^2 f”(x) + o(h^2) $$
これを積分に代入し、カーネルの3条件（$\int K=1,\ \int uK=0,\ \int u^2 K=\mu_2(K)$）を順に当てはめると、各項がきれいに片付きます。第1項は $f(x)\int K(u)du = f(x)$、第2項は $hf'(x)\int uK(u)du = 0$、第3項は $\frac{1}{2}h^2 f”(x)\int u^2K(u)du = \frac{1}{2}h^2\mu_2(K)f”(x)$ です。
$$ \mathbb{E}[\hat{f}_h(x)] = f(x) + \frac{1}{2}h^2 \mu_2(K)\, f”(x) + o(h^2) $$
したがってバイアスは、期待値から真値を引いて、
$$ \begin{equation} \mathrm{Bias}[\hat{f}_h(x)] = \mathbb{E}[\hat{f}_h(x)] – f(x) = \frac{1}{2}h^2 \mu_2(K)\, f”(x) + o(h^2) \end{equation} $$
と求まります。ここから2つのことが読み取れます。バイアスは $h^2$ に比例して大きくなる（$h$ を広げるほどバイアス増）こと、そして密度の曲率 $f”(x)$ に比例すること。山や谷の頂点（$f”$ が大きい場所）ほど大きく平滑化されて潰れる、という直感とぴったり一致します。
この図はバイアスを目に見える形にしたものです。破線が真の密度、青線がKDEの期待値 $\mathbb{E}[\hat{f}_h]$（=真の密度とカーネルの畳み込み）です。山の頂点（$f”<0$）ではオレンジのぶんだけ低く潰れ、谷の底（$f''>0$）では緑のぶんだけ持ち上がっています。ずれの向きと大きさが、ちょうど曲率 $f”(x)$ に従っている——式 $\frac{1}{2}h^2\mu_2 f”$ のとおりです。
分散の導出
次に分散です。独立な $n$ 項の平均の分散は1項の分散の $1/n$ なので、
$$ \mathrm{Var}[\hat{f}_h(x)] = \frac{1}{n}\,\mathrm{Var}\!\left[\frac{1}{h}K\!\left(\frac{x-X}{h}\right)\right] = \frac{1}{n}\left( \frac{1}{h^2}\mathbb{E}\!\left[K^2\!\left(\tfrac{x-X}{h}\right)\right] – \frac{1}{h^2}\Big(\mathbb{E}\big[K(\tfrac{x-X}{h})\big]\Big)^2 \right) $$
第1項を、バイアスのときと同じ置換 $u=(t-x)/h$ で書き直します。$\frac{1}{h^2}\mathbb{E}[K^2] = \frac{1}{h}\int K^2(u)f(x+hu)\,du$ となり、$h\to 0$ では $f(x+hu)\approx f(x)$ なので、
$$ \frac{1}{h^2}\mathbb{E}\!\left[K^2\!\left(\tfrac{x-X}{h}\right)\right] \approx \frac{1}{h}\, f(x) \int K^2(u)\,du = \frac{f(x)\, R(K)}{h} $$
と近似できます。ここで $R(K) = \int K^2(u)\,du$ と置きました。一方、第2項の $\frac{1}{h^2}(\mathbb{E}[K])^2$ は $\mathbb{E}[K(\cdot)/h]\to f(x)$ より $O(1)$ の有限量で、$1/h$ で発散する第1項に比べれば無視できます。主要項だけ残すと、
$$ \begin{equation} \mathrm{Var}[\hat{f}_h(x)] = \frac{f(x)\, R(K)}{nh} + o\!\left(\frac{1}{nh}\right) \end{equation} $$
を得ます。分散は $1/(nh)$ に比例します。$h$ を小さくすると分散が増え、$n$ を増やすかバンド幅を広げると分散が減る——「狭い山ほどデータのばらつきに敏感」という直感どおりです。
分散とは「データを取り直すたびに推定がどれだけ暴れるか」です。同じ真の密度から200点を30回取り直し、それぞれのKDEを薄い線で重ねました。左（$h=0.15$）は線がばらばらに散らばっています——分散が大きい。一方で平均的には真の密度をよく追えています（バイアスは小さい）。右（$h=0.8$）は線がぴたりと束になり安定（分散が小さい）ですが、束ごと真の密度から内側にずれています（バイアスが大きい）。左右でバイアスと分散がちょうど逆——これがトレードオフの正体です。
ここで決定的な綱引きが見えます。バイアスは $h^2$ で増え、分散は $1/(nh)$ で増える。$h$ を大きくすればバイアスが、小さくすれば分散が悪化する。両方を同時に小さくはできません。この綱引きをひとつの指標にまとめたのが、次の AMISE です。
AMISEと最適バンド幅
各点での誤差を全域で集計するため、平均積分二乗誤差（MISE） を考えます。これは「真値とのずれの二乗」を期待値・全域積分したもので、バイアスの二乗と分散の積分に分解できます（二乗誤差＝バイアス²＋分散、という標準的な分解です）。
$$ \mathrm{MISE}(h) = \int \mathbb{E}\big[(\hat{f}_h(x)-f(x))^2\big]\,dx = \int \mathrm{Bias}^2\,dx + \int \mathrm{Var}\,dx $$
前の節で求めたバイアスと分散の主要項を代入します。バイアス二乗の積分は $\int \big(\frac{1}{2}h^2\mu_2 f”(x)\big)^2 dx = \frac{1}{4}h^4 \mu_2(K)^2 \int f”(x)^2 dx$、分散の積分は $\int \frac{f(x)R(K)}{nh}dx = \frac{R(K)}{nh}$（$\int f\,dx=1$ より）です。これらをまとめた漸近形を AMISE（漸近 MISE） と呼びます。
$$ \begin{equation} \mathrm{AMISE}(h) = \frac{1}{4}h^4 \mu_2(K)^2\, R(f”) + \frac{R(K)}{nh} \end{equation} $$
ここで $R(f”) = \int f”(x)^2 dx$ は真の密度の「ごつごつ具合」を表す量です。第1項（バイアス由来）は $h^4$ で増え、第2項（分散由来）は $1/(nh)$ で減る。この和を最小にする $h$ が最適バンド幅です。$h$ で微分してゼロと置きます。
$$ \frac{d\,\mathrm{AMISE}}{dh} = h^3 \mu_2(K)^2 R(f”) – \frac{R(K)}{nh^2} = 0 $$
$h$ について解くために両辺に $nh^2$ を掛けて整理すると $h^5 = \dfrac{R(K)}{n\,\mu_2(K)^2 R(f”)}$ となり、
$$ \begin{equation} h_{\mathrm{opt}} = \left(\frac{R(K)}{n\,\mu_2(K)^2\, R(f”)}\right)^{1/5} \propto n^{-1/5} \end{equation} $$
が得られます。ここから2つの大事な事実が読めます。第一に、最適バンド幅はサンプル数とともに $n^{-1/5}$ でゆっくり縮む——データが増えるほど山を細くしてよい。第二に、これを AMISE に戻すと最小誤差は $O(n^{-4/5})$ となり、パラメトリック推定の $O(n^{-1})$ よりわずかに遅い。これがノンパラメトリックの「自由さの代償」です。
この $n^{-1/5}$ という法則を実験で確かめたのが上の図です。サンプル数 $n$ を $50$ から $2000$ まで変え、各 $n$ で誤差を最小にするバンド幅を実験的に求めてプロットしました（青点）。両対数グラフ上で、点が理論の傾き $-1/5$ の直線（オレンジ破線）にきれいに乗っています。データを10倍にしても最適バンド幅は $10^{-1/5}\approx 0.63$ 倍にしか縮まない——「もっとデータがあれば一気に細かく見られる」とはいかない、というノンパラメトリックの宿命がはっきり見えます。
この「誤差が $h$ に対してどう動くか」は、後の実装で実際に図にして確かめます（下の MISE 曲線の節）。理論上は、誤差を $h$ の関数として描くと、最適バンド幅のところで底になる谷型のカーブになるはずです。
ところが $h_{\mathrm{opt}}$ は真の密度の曲率 $R(f”)$ を含んでいて、それが分かっていれば苦労はありません（鶏と卵の問題）。そこで実務では、$f$ を正規分布と仮定して $R(f”)$ を見積もる近似を使います。次の Silverman のルールがそれです。
Silvermanの目安公式
真の密度をガウスカーネルのもとで正規分布 $N(\mu,\sigma^2)$ と仮定して $R(f”)$ を計算し、$h_{\mathrm{opt}}$ に代入すると、有名な Silverman の目安公式（rule of thumb） が導かれます。
$$ \begin{equation} h_{\mathrm{Silverman}} = 0.9\,\min\!\left(\hat{\sigma},\ \frac{\mathrm{IQR}}{1.34}\right) n^{-1/5} \end{equation} $$
$\hat{\sigma}$ は標本標準偏差、$\mathrm{IQR}$ は四分位範囲です。標準偏差の代わりに $\mathrm{IQR}/1.34$ との最小を取るのは、外れ値で $\hat{\sigma}$ が膨らんでバンド幅が過大になるのを防ぐためです。係数 $0.9$ は、単峰の正規分布で導いた $1.06$ を、多峰分布でも平滑化しすぎないよう少し小さくした実用値です。万能ではありませんが、出発点としてとても便利です。理論が出そろったので、Python で動かして確かめましょう。
Pythonでの実装
準備とKDEのスクラッチ実装
まず、定義式 $\hat{f}_h(x)=\frac{1}{nh}\sum_i K((x-x_i)/h)$ をそのまま実装します。カーネルはガウスを使います。真の密度として2つの山を持つ混合正規分布を用意し、そこから標本を取って KDE で復元してみます。図の軸ラベルを日本語にするため、最初にフォントだけ設定しておきます。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# 日本語ラベルを表示するためのフォント設定（環境にあるものを自動選択）
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
rng = np.random.default_rng(0)
# 真の密度：2つの山を持つ混合正規分布
def true_pdf(x):
g = lambda x, m, s: np.exp(-0.5*((x-m)/s)**2) / (s*np.sqrt(2*np.pi))
return 0.6*g(x, -2.0, 0.7) + 0.4*g(x, 2.0, 1.0)
# 混合分布からのサンプリング（n=200）
n = 200
comp = rng.random(n) < 0.6
data = np.where(comp, rng.normal(-2.0, 0.7, n), rng.normal(2.0, 1.0, n))
def gaussian_kernel(u):
return np.exp(-0.5*u**2) / np.sqrt(2*np.pi)
def kde(x_grid, data, h):
# x_grid:(m,), data:(n,) -> 各グリッド点で密度を返す
u = (x_grid[:, None] - data[None, :]) / h # (m, n)
return gaussian_kernel(u).sum(axis=1) / (len(data) * h)
xs = np.linspace(-6, 6, 400)
fhat = kde(xs, data, h=0.5)
plt.figure(figsize=(8, 4.6))
plt.plot(xs, true_pdf(xs), "--", color="gray", lw=2, label="真の密度（正解）")
plt.plot(xs, fhat, color="#2c7fb8", lw=2.5, label="KDE推定（h=0.5）")
plt.plot(data, np.zeros_like(data)-0.006, "|", color="gray", alpha=0.6, label="観測データ200点")
plt.xlabel("データの値"); plt.ylabel("密度"); plt.legend()
plt.title("200個のデータだけから密度を復元する")
plt.tight_layout(); plt.show()
この図から、わずか200点の標本で2つの山がきちんと復元できているのがわかります。各データ点（下の縦棒）の上にガウスの山が積み上がり、点が密集する $x\approx-2$ と $x\approx 2$ で密度が高くなっています。真の密度（破線）と KDE（実線）はよく重なっていて、ヒストグラムのような階段状の段差もありません。次は、この結果を左右する最重要パラメータ——バンド幅 $h$ の効果を見ます。
バンド幅による過小平滑化と過剰平滑化
理論で予測した「$h$ 小→分散大（ギザギザ）、$h$ 大→バイアス大（のっぺり）」を実際に確かめます。3つのバンド幅で KDE を描き比べましょう。
specs = [(0.1, "狭すぎ：ギザギザ（過小平滑化）", "#e6550d"),
(0.5, "ちょうど良い", "#2ca25f"),
(2.0, "広すぎ：山が潰れる（過剰平滑化）", "#2c7fb8")]
fig, axes = plt.subplots(1, 3, figsize=(15, 4.4), sharey=True)
for ax, (h, t, c) in zip(axes, specs):
ax.plot(xs, true_pdf(xs), "--", color="gray", lw=2, label="真の密度")
ax.plot(xs, kde(xs, data, h), color=c, lw=2.5, label=f"KDE（h={h}）")
ax.plot(data, np.zeros_like(data)-0.006, "|", color="gray", alpha=0.4)
ax.set_title(t); ax.set_xlabel("データの値"); ax.legend()
axes[0].set_ylabel("密度")
plt.tight_layout(); plt.show()
3枚の比較が、バイアスと分散のトレードオフを目で見せてくれます。左（$h=0.1$）は山が細すぎて、推定がデータの偶然のばらつきを拾い、いくつもの偽の小さな山が出ています——過小平滑化（分散大）です。右（$h=2.0$）は山が広すぎて2つの峰が溶け合い、1つのなだらかな丘に潰れています——過剰平滑化（バイアス大）で、本来の2峰構造が消えています。中央（$h=0.5$）が両者のバランスが取れた良い推定です。では、この「ちょうどよい $h$」は理論どおり谷型の誤差カーブの底にあるのでしょうか。次に確かめます。
バンド幅と誤差(MISE)の関係を実験で確かめる
AMISE の節で「誤差を $h$ の関数として描くと谷型になり、その底が最適バンド幅」と予言しました。これを実験で再現します。いろいろな $h$ について、KDE と真の密度の二乗誤差を全域で積分した量（ISE）を、データを取り直しながら何度も測って平均します。これが MISE の近似です。
hs = np.logspace(np.log10(0.05), np.log10(3.0), 40)
n_trials = 30
mise = np.zeros_like(hs)
ft = true_pdf(xs)
for _ in range(n_trials):
cc = rng.random(n) < 0.6
d = np.where(cc, rng.normal(-2.0, 0.7, n), rng.normal(2.0, 1.0, n))
for j, h in enumerate(hs):
mise[j] += np.trapz((kde(xs, d, h) - ft)**2, xs) # ISEを積分
mise /= n_trials
h_best = hs[mise.argmin()]
print(f"実験的に最適なバンド幅: {h_best:.3f}")
plt.figure(figsize=(8, 4.8))
plt.plot(hs, mise, "o-", color="#2c7fb8", ms=4)
plt.axvline(h_best, color="#2ca25f", ls="--", label=f"最適バンド幅 ≈ {h_best:.2f}")
plt.xscale("log")
plt.xlabel("バンド幅 h（対数軸）"); plt.ylabel("平均積分二乗誤差 MISE")
plt.title("バンド幅と推定誤差(MISE)：谷の底が最適"); plt.legend()
plt.tight_layout(); plt.show()
この谷型のカーブこそ、理論で導いた AMISE の姿そのものです。左側（$h$ が小さい領域）は分散が支配して誤差が大きく、右側（$h$ が大きい領域）はバイアスが支配して誤差が大きい。その間の $h\approx 0.33$ で誤差が最小になっています。式の上だけだった「バイアスと分散の綱引き」が、1本のなめらかな谷として目に見えました。次は、この最適値を真の密度を知らずに当てにいく Silverman の公式を試します。
Silvermanの公式による自動バンド幅選択
導出した Silverman の目安公式を実装し、それが返すバンド幅を確かめます。
def silverman_bandwidth(data):
n = len(data)
sigma = np.std(data, ddof=1)
iqr = np.subtract(*np.percentile(data, [75, 25]))
spread = min(sigma, iqr / 1.34)
return 0.9 * spread * n ** (-1/5)
h_sil = silverman_bandwidth(data)
print(f"Silvermanのバンド幅: {h_sil:.3f}")
実行すると、Silverman のバンド幅は $h\approx 0.69$ を返します。先ほど実験で見つけた最適値 $0.33$ よりやや広めです。これは偶然ではありません。この公式は内部で「真の密度は単峰の正規分布」と仮定しているため、今回のように峰が2つある分布では、本来より広め＝過剰平滑化気味に出る傾向があります。とはいえデータの標準偏差とサンプル数だけから、追加の調整なしに同じオーダーの妥当な値が一発で得られるのは強力です。多峰性が疑われるときは交差検証で $h$ を選ぶのが安全、と覚えておきましょう。次は、KDE を2次元へ広げます。
2次元への拡張
KDE は多次元にも自然に広がります。各点に2次元のガウスの山を置き、足し合わせるだけです。ここでは2つのクラスタを持つ2次元データで密度を推定し、等高線で可視化します。
# 2つのクラスタを持つ2次元データ
n2 = 300
c1 = rng.normal([-1.5, -1.0], 0.6, size=(n2//2, 2))
c2 = rng.normal([1.5, 1.0], 0.8, size=(n2//2, 2))
pts = np.vstack([c1, c2])
def kde_2d(grid, data, h):
diff = grid[:, None, :] - data[None, :, :] # (G, n, 2)
sq = (diff**2).sum(axis=2) / h**2 # (G, n)
k = np.exp(-0.5*sq) / (2*np.pi) # 2次元標準ガウス
return k.sum(axis=1) / (len(data) * h**2)
gx, gy = np.meshgrid(np.linspace(-4, 4, 120), np.linspace(-4, 4, 120))
grid = np.column_stack([gx.ravel(), gy.ravel()])
dens = kde_2d(grid, pts, h=0.5).reshape(gx.shape)
plt.figure(figsize=(6.8, 5.6))
cf = plt.contourf(gx, gy, dens, levels=20, cmap="viridis")
plt.scatter(pts[:, 0], pts[:, 1], s=8, c="white", alpha=0.45)
plt.colorbar(cf, label="推定された密度")
plt.xlabel("特徴 1"); plt.ylabel("特徴 2"); plt.title("2次元のカーネル密度推定")
plt.tight_layout(); plt.show()
等高線図から、2つのクラスタに対応する2つの密度のピークがはっきり浮かび上がっています。白い点（データ）が密集する2か所で密度が最大になり、その間で密度が落ちています。1次元のときと同じ「各点に山を置いて重ねる」原理が、次元を上げてもそのまま通用するわけです。ただし高次元では、空間が広大になって点がスカスカになる「次元の呪い」のため、バンド幅選びが難しくなる点には注意が必要です。最後に、自作した KDE がライブラリと一致するかを確かめます。
scipyとの比較で実装を検証
自作の KDE が正しいかを、scipy.stats.gaussian_kde と比べて確かめます。scipy はバンド幅をデータの標準偏差を単位として指定するので、自作の $h$ を $\hat{\sigma}$ で割って渡します。
from scipy.stats import gaussian_kde
kde_scipy = gaussian_kde(data, bw_method=h_sil / np.std(data, ddof=1))
fhat_scipy = kde_scipy(xs)
fhat_mine = kde(xs, data, h_sil)
print("最大絶対誤差:", np.max(np.abs(fhat_mine - fhat_scipy)))
plt.figure(figsize=(8, 4.6))
plt.plot(xs, fhat_mine, color="#2c7fb8", lw=4, label="自作KDE")
plt.plot(xs, fhat_scipy, "--", color="#e6550d", lw=2, label="scipy.stats.gaussian_kde")
plt.xlabel("データの値"); plt.ylabel("密度"); plt.legend()
plt.title("自作KDEとscipyが一致することの確認")
plt.tight_layout(); plt.show()
2本の曲線はほぼ完全に重なり、最大絶対誤差も $10^{-16}$ 程度のごく小さな値にとどまります。これは自作のスクラッチ実装が、定義に忠実で正しく動いている証拠です。理論どおりの推定が、わずか数行のコードで実現できることが確かめられました。
まとめ
本記事では、カーネル密度推定を理論から導出し、図とともに Python で実装しました。
- 動機: ヒストグラムのビン境界の恣意性と階段状の段差を解消するため、各データ点に山（カーネル）を置いて重ね合わせる。$\hat{f}_h(x)=\frac{1}{nh}\sum_i K((x-x_i)/h)$。
- バイアス: $\frac{1}{2}h^2\mu_2(K)f”(x)$。$h^2$ と密度の曲率に比例し、山や谷を潰す。
- 分散: $\frac{f(x)R(K)}{nh}$。$1/(nh)$ に比例し、$h$ が小さいほど大きい。
- トレードオフと AMISE: バイアス（$h^4$）と分散（$1/nh$）の和を最小にする最適バンド幅は $h_{\mathrm{opt}}\propto n^{-1/5}$。実験でも誤差は谷型になり、底が最適値に一致した。
- Silverman の公式: $0.9\min(\hat\sigma, \mathrm{IQR}/1.34)\,n^{-1/5}$。正規分布仮定の実用的なバンド幅の出発点。多峰分布ではやや広めに出る。
- 実装: スクラッチ KDE でバンド幅の効果、MISE の谷、2次元 KDE、scipy との一致を確かめた。
カーネル密度推定は「分布の形を仮定せず、データに語らせる」ノンパラメトリック統計の入り口であり、異常検知・可視化・ベイズ推定の土台になります。バンド幅選びがバイアスと分散のトレードオフそのものだ、という視点を持っておくと、機械学習の正則化やモデル選択とも一本の線でつながって見えてきます。
次のステップとして、以下の記事も参考にしてください。
- 混合分布とは？混合正規分布の理論とEM推定をわかりやすく解説 — KDE をパラメトリックに置き換えた混合モデル
- バイアス・バリアンスのトレードオフ（統計） — バンド幅選びの本質
- 正規分布とは？定義やグラフの見方を徹底解説 — ガウスカーネルの基礎
参考文献
- B. W. Silverman, “Density Estimation for Statistics and Data Analysis,” Chapman & Hall, 1986.
- M. P. Wand and M. C. Jones, “Kernel Smoothing,” Chapman & Hall, 1995.
- D. W. Scott, “Multivariate Density Estimation: Theory, Practice, and Visualization,” Wiley, 1992.