第05回カーネル密度推定（KDE）¶

科目： ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕
担当： 小野原彩香
対象： 3年生（後期）

学習目標¶

カーネル密度推定（KDE）の数学的な原理を説明できる
バンド幅（h）がKDEの結果に与える影響を理解し、適切に選択できる
gdal_grid コマンドを使ってポイントデータからKDEラスターを生成できる
scikit-learnの KernelDensity を用いてPythonでKDEを計算できる
foliumまたはmatplotlibでヒートマップを可視化・出力できる
発展：KDEの式が確率密度になる理由と、バンド幅 \(h\) の役割を、自分で導ける

1. 密度推定とは¶

1.1 問題の設定¶

地理空間データを扱う際、個々のポイント（例：LLMへのアクセスログの位置情報）がどの地域に集中しているかを把握したいことがよくあります。しかし、ポイントの数が多くなると、単純な散布図では密度の高低が見えにくくなります。

密度推定（Density Estimation）とは、観測されたポイントデータから背後にある確率密度関数を推定する手法です。密度推定を使うと、「この地域にはどのくらいの頻度でアクセスが集中しているか」を連続的な曲面として表現できます。

1.2 ヒストグラムとの違い¶

手法	特徴
ヒストグラム	セル（ビン）の境界に依存、不連続
KDE	滑らかな連続曲面、境界の影響を受けない

2. カーネル密度推定（KDE）の原理¶

2.1 KDEの数式¶

\(n\) 個の観測点 \(x_1, x_2, \ldots, x_n\) が与えられたとき、点 \(x\) における密度推定値は次の式で定義されます。

\[\hat{f}(x) = \frac{1}{nh} \sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right)\]

ここで： - \(\hat{f}(x)\)：点 \(x\) における推定密度 - \(n\)：観測点の数 - \(h\)：バンド幅（bandwidth）。平滑化の程度を制御するパラメータ - \(K(\cdot)\)：カーネル関数。各観測点の影響の形を定義する

導出①：なぜ \(\dfrac{1}{nh}\) が前に付くのか ―― \(\hat{f}\) が確率密度になる理由

KDEの式は「観測点ごとに小さな山を置いて、平均する」という操作です。前に付く \(\frac{1}{nh}\) は、その結果がちゃんと確率密度になるための係数です。

観測点 \(x_i\) が1つ作る山だけを取り出します。

\[g_i(x) = \frac{1}{h}K\!\left(\frac{x-x_i}{h}\right)\]

この面積を計算します。\(u = \dfrac{x-x_i}{h}\) と置くと \(dx = h\,du\) なので、

\[\int_{-\infty}^{\infty} g_i(x)\,dx = \int_{-\infty}^{\infty} \frac{1}{h}K(u)\cdot h\,du = \int_{-\infty}^{\infty} K(u)\,du = 1\]

\(h\) がきれいに約分されて消えます。カーネル \(K\) 自体の面積が1だから（ガウシアンカーネルは標準正規分布そのものです）、山1つの面積も \(h\) によらず必ず1になります。

\(\hat{f}\) はこの山を \(n\) 個ぶん足して \(n\) で割ったもの、つまり平均です。

\[\hat{f}(x) = \frac{1}{n}\sum_{i=1}^{n} g_i(x) = \frac{1}{nh}\sum_{i=1}^{n} K\!\left(\frac{x-x_i}{h}\right)\]

面積1のものを \(n\) 個平均したのだから、全体の面積も1になります。これで \(\hat{f}\) は確率密度の条件（非負・面積1）を満たします。

数値で確かめる：観測点を \(\{1.0,\ 2.0,\ 2.5,\ 4.0,\ 7.0\}\) の5点として、ガウシアンカーネルで数値積分すると \(\int K(u)\,du = 1.0000000000\)、\(\int \hat{f}(x)\,dx\) は \(h = 0.3,\ 1.0,\ 3.0\) のいずれでも \(1.00000000\) になります。

使った数学：置換積分（微I）。平均が線形であること。

2.2 ガウシアンカーネル¶

最もよく使われるカーネル関数はガウシアン（正規分布）カーネルです。

\[K(u) = \frac{1}{\sqrt{2\pi}} \exp\!\left(-\frac{u^2}{2}\right)\]

2次元（緯度・経度）の場合は、各次元にガウシアンカーネルを適用します。

\[\hat{f}(x, y) = \frac{1}{n h_x h_y} \sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h_x}\right) K\!\left(\frac{y - y_i}{h_y}\right)\]

2.3 バンド幅（h）の役割¶

バンド幅 \(h\) はKDEで最も重要なパラメータです。

h が小さい場合：各観測点の影響が局所的になり、細かいピークが多く現れる（過適合）
h が大きい場合：観測点の影響が広く広がり、全体的に平滑化される（過平滑化）

適切な \(h\) を選ぶことが、KDEの品質を大きく左右します。

導出②：面積は1で固定されている ―― だから \(h\) は「高さと広がりの配分」を決める

導出①で、\(h\) をどう変えても \(\hat{f}\) の面積は1のままだと分かりました。ここから、バンド幅のトレードオフが式の帰結として出てきます。

面積が変えられない以上、山を細くすれば高くなるしかなく、広げれば低くなるしかない。\(h\) が決めているのは総量ではなく、その配分です。

\(h\) を小さくする → 各観測点の山が細く高くなる → 観測点の位置に鋭いピークが立つ（過適合）
\(h\) を大きくする → 各観測点の山が広く低くなる → 隣どうしが重なって溶け合う（過平滑化）

数値で確かめる：同じ5点 \(\{1.0,\ 2.0,\ 2.5,\ 4.0,\ 7.0\}\) で \(h\) だけを変えると、

\(h\)	\(\int \hat{f}(x)\,dx\)	\(\hat{f}\) の最大値
0.3	1.00000000	0.3759
1.0	1.00000000	0.2097
3.0	1.00000000	0.1089

面積は動かず、高さだけが動いています。

カーネル \(K\) をガウシアン以外（エパネチニコフなど）に変えても、どれも面積1で形が少し違うだけなので、結果は \(h\) を変えたときほどには変わりません。これはよく知られた事実ですが、どのカーネルがどれだけ効率的かを比べる理論は、この授業で扱う範囲外です。

使った数学：導出①の面積の保存。

2.4 Silvermanの経験則（Rule of Thumb）¶

データ数 \(n\)、標準偏差 \(\sigma\) のとき、Silvermanの経験則によるバンド幅の推定式は次のとおりです。

\[h = 0.9 \cdot \min\!\left(\sigma, \frac{\mathrm{IQR}}{1.34}\right) \cdot n^{-1/5}\]

ここで IQR は四分位範囲（Interquartile Range）です。この式は自動的に適切なバンド幅を計算するため、初めてKDEを試みる際の出発点として広く使われています。

導出③：Silvermanの式に \(n^{-1/5}\) が出てくる意味

この式は2つの部分でできています。

\[h = \underbrace{0.9 \cdot \min\!\left(\sigma, \frac{\mathrm{IQR}}{1.34}\right)}_{\text{データの広がり}} \cdot \underbrace{n^{-1/5}}_{\text{データ数}}\]

前半（広がり）：バンド幅はデータのスケールに比例すべきです。メートルで測ったデータとキロメートルで測ったデータでは、適切な \(h\) も1000倍違うからです。\(\sigma\) と \(\mathrm{IQR}/1.34\) の小さいほうを採るのは、外れ値に引きずられないためです。正規分布では \(\mathrm{IQR}\) が \(\sigma\) の約1.34倍になるので、1.34で割ると両者が同じ尺度に揃います。外れ値があると \(\sigma\) だけが大きくなるため、小さいほうを採れば影響を避けられます。

後半（データ数）：\(h\) を小さくすると細かい形は追えますが、山が細いぶん少数の点に左右されて推定が不安定になります。逆に \(h\) を大きくすると安定しますが、形がなまります。この2つは \(h\) に対して逆向きに動くので、釣り合う点があります。その釣り合いを取ると \(h\) は \(n^{-1/5}\) に比例する、というのが指数の出どころです。

指数が \(-1/5\) と小さいので、データを増やしてもバンド幅はあまり小さくなりません。

データ数の変化	\(h\) の倍率
10倍	0.6310
100倍	0.3981

データを10倍集めても \(h\) は約63%にしかなりません。「データを増やせば細かい構造が見える」は、思ったほど速くは効かないということです。

ずれとばらつきを具体的に式に書いて釣り合いを解く計算（平均二乗積分誤差の漸近展開）は、この授業で扱う範囲外です。ここでは「2つが逆向きに動くから最適値がある」「その結果が \(n^{-1/5}\) になる」までを押さえてください。

使った数学：比例関係と指数の計算。釣り合いの考え方。

3. gdalによるKDE¶

3.1 gdal_grid の概要¶

gdal_grid はポイントデータを規則的なグリッド（ラスター）に変換するGDALのコマンドラインツールです。補間アルゴリズムの一つとして、ガウシアンカーネルを用いた移動平均（moving_average）や逆距離加重法（invdist）なども利用可能です。

3.2 入力データの準備¶

CSVファイルをVRTに変換する方法が最も簡単です。

# ポイントデータ（CSV）の例
# llm_access_log_points.csv
# lon,lat,count
# 141.35,43.07,5
# 141.36,43.08,12
# ...

# CSVをGeoJSONに変換（OGRを使用）
ogr2ogr \
  -f GeoJSON \
  llm_access_points.geojson \
  llm_access_log_points.csv \
  -oo X_POSSIBLE_NAMES=lon \
  -oo Y_POSSIBLE_NAMES=lat \
  -a_srs EPSG:4326

3.3 gdal_grid によるKDEラスター生成¶

# ガウシアン移動平均によるKDE的ラスター生成
# radius1, radius2: X/Y方向の検索半径（度単位）
# angle: 回転角（0で軸平行）
# min_points: 最低点数（0で常に計算）
# nodata: NoData値

gdal_grid \
  -zfield "count" \
  -a moving_average:radius1=0.05:radius2=0.05:angle=0:min_points=0:nodata=-9999 \
  -txe 140.9 141.8 \
  -tye 42.7 43.4 \
  -outsize 256 256 \
  -of GTiff \
  -ot Float32 \
  llm_access_points.geojson \
  kde_output.tif

# 結果を確認
gdalinfo kde_output.tif

3.4 ラスターの投影変換と色付け¶

# 地図表示用にカラーランプを適用
gdaldem color-relief kde_output.tif color_ramp.txt kde_colored.tif

# color_ramp.txt の例
# 値  R    G    B
# 0   255  255  255
# 5   255  255  0
# 20  255  128  0
# 50  255  0    0
# 100 128  0    0

4. PythonによるKDE（scipy/scikit-learn）¶

4.1 scikit-learnのKernelDensity¶

import numpy as np
import pandas as pd
from sklearn.neighbors import KernelDensity
import matplotlib.pyplot as plt

# ログデータの読み込み
df = pd.read_csv("llm_access_log_points.csv")

# 座標を配列に変換（緯度・経度）
coords = df[["lon", "lat"]].values  # shape: (n, 2)

# Silvermanの経験則でバンド幅を推定
n = len(coords)
sigma_lon = np.std(coords[:, 0])
sigma_lat = np.std(coords[:, 1])
sigma = np.mean([sigma_lon, sigma_lat])
iqr = np.percentile(coords, 75, axis=0) - np.percentile(coords, 25, axis=0)
iqr_mean = np.mean(iqr) / 1.34

h_silverman = 0.9 * min(sigma, iqr_mean) * n ** (-0.2)
print(f"Silvermanのバンド幅推定値: {h_silverman:.5f} 度")

# KDEモデルの構築
kde = KernelDensity(kernel="gaussian", bandwidth=h_silverman)
kde.fit(coords)

# 評価グリッドの作成
lon_min, lon_max = 140.9, 141.8
lat_min, lat_max = 42.7, 43.4
grid_size = 200

lon_grid = np.linspace(lon_min, lon_max, grid_size)
lat_grid = np.linspace(lat_min, lat_max, grid_size)
lon_mesh, lat_mesh = np.meshgrid(lon_grid, lat_grid)
grid_points = np.column_stack([lon_mesh.ravel(), lat_mesh.ravel()])

# 密度の計算（log スケールで返ってくるので exp で戻す）
log_density = kde.score_samples(grid_points)
density = np.exp(log_density).reshape(grid_size, grid_size)

print(f"密度の最大値: {density.max():.6f}")
print(f"密度の最小値: {density.min():.6f}")

4.2 scipyを使った方法¶

from scipy.stats import gaussian_kde

# scipy版（1次元・多次元対応）
coords_T = coords.T  # shape: (2, n)に転置が必要
kde_scipy = gaussian_kde(coords_T, bw_method="silverman")

# グリッド上で評価
density_scipy = kde_scipy(grid_points.T).reshape(grid_size, grid_size)

print("scipy版 KDE 計算完了")

4.3 GeoTIFFへの書き出し¶

import rasterio
from rasterio.transform import from_bounds
from rasterio.crs import CRS

# GeoTIFFとして保存
transform = from_bounds(
    lon_min, lat_min, lon_max, lat_max,
    grid_size, grid_size
)

with rasterio.open(
    "kde_python_output.tif",
    "w",
    driver="GTiff",
    height=grid_size,
    width=grid_size,
    count=1,
    dtype=density.dtype,
    crs=CRS.from_epsg(4326),
    transform=transform,
) as dst:
    # ラスターは上下反転して保存（GeoTIFF の Y 軸は上が大きい）
    dst.write(np.flipud(density), 1)

print("GeoTIFF を kde_python_output.tif に保存しました")

5. ヒートマップの可視化¶

5.1 matplotlibによる可視化¶

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy.crs as ccrs  # インストール必要: pip install cartopy

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# --- バンド幅の比較 ---
for ax, bw, title in zip(
    axes,
    [h_silverman * 0.3, h_silverman],
    [f"小さいバンド幅 (h={h_silverman*0.3:.4f})", f"Silverman則 (h={h_silverman:.4f})"]
):
    kde_temp = KernelDensity(kernel="gaussian", bandwidth=bw)
    kde_temp.fit(coords)
    log_d = kde_temp.score_samples(grid_points)
    d = np.exp(log_d).reshape(grid_size, grid_size)

    im = ax.imshow(
        np.flipud(d),
        extent=[lon_min, lon_max, lat_min, lat_max],
        aspect="auto",
        cmap="hot_r",
        origin="upper"
    )
    ax.scatter(coords[:, 0], coords[:, 1], s=2, c="blue", alpha=0.3, label="観測点")
    ax.set_title(title)
    ax.set_xlabel("経度")
    ax.set_ylabel("緯度")
    ax.legend(fontsize=8)
    plt.colorbar(im, ax=ax, label="密度")

plt.suptitle("KDE バンド幅の比較（LLMアクセスログ）", fontsize=14)
plt.tight_layout()
plt.savefig("kde_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
print("PNG を kde_comparison.png に保存しました")

5.2 foliumによるインタラクティブヒートマップ¶

import folium
from folium.plugins import HeatMap

# foliumの HeatMap プラグインを使用
m = folium.Map(
    location=[43.05, 141.35],  # 中心座標（札幌付近）
    zoom_start=11,
    tiles="CartoDB positron"
)

# HeatMap には [lat, lon, weight] のリストを渡す
heat_data = [
    [row["lat"], row["lon"], row.get("count", 1)]
    for _, row in df.iterrows()
]

HeatMap(
    heat_data,
    radius=20,         # ピクセル半径
    blur=15,           # ブラー強度
    max_zoom=14,
    gradient={
        "0.2": "blue",
        "0.5": "lime",
        "0.8": "orange",
        "1.0": "red"
    }
).add_to(m)

# 凡例を追加
folium.map.Marker(
    [43.35, 141.7],
    icon=folium.DivIcon(
        html='<div style="font-size:12px;background:white;padding:5px;border:1px solid gray;">'
             'LLMアクセス密度ヒートマップ<br>（カーネル密度推定）</div>'
    )
).add_to(m)

m.save("kde_heatmap.html")
print("インタラクティブマップを kde_heatmap.html に保存しました")

6. ハンズオン手順¶

事前準備¶

# 必要ライブラリのインストール
pip install scikit-learn scipy rasterio folium pandas numpy matplotlib

# GDALのインストール（macOS）
brew install gdal

# データの確認
ls -lh data/
# 期待するファイル: llm_access_log_points.csv

ステップ1: データの確認¶

import pandas as pd

df = pd.read_csv("data/llm_access_log_points.csv")
print(df.head())
print(f"\n総レコード数: {len(df)}")
print(f"緯度範囲: {df['lat'].min():.4f} ～ {df['lat'].max():.4f}")
print(f"経度範囲: {df['lon'].min():.4f} ～ {df['lon'].max():.4f}")
print(f"\n欠損値:\n{df.isnull().sum()}")

ステップ2: KDEの実行とGeoTIFF出力¶

上記セクション4.1〜4.3のコードを順番に実行してください。

ステップ3: PNG出力¶

上記セクション5.1のコードを実行してください。

ステップ4: foliumマップ出力¶

上記セクション5.2のコードを実行して、ブラウザで kde_heatmap.html を開いて確認してください。

ステップ5: バンド幅の変更実験¶

# バンド幅を変えて結果を比較する
bandwidths = [0.005, 0.01, 0.05, 0.1]

fig, axes = plt.subplots(1, len(bandwidths), figsize=(20, 5))
for ax, bw in zip(axes, bandwidths):
    kde_test = KernelDensity(kernel="gaussian", bandwidth=bw)
    kde_test.fit(coords)
    log_d = kde_test.score_samples(grid_points)
    d = np.exp(log_d).reshape(grid_size, grid_size)
    ax.imshow(np.flipud(d), extent=[lon_min, lon_max, lat_min, lat_max],
              aspect="auto", cmap="YlOrRd", origin="upper")
    ax.set_title(f"h = {bw}")
    ax.axis("off")

plt.suptitle("バンド幅の違いによるKDE比較")
plt.tight_layout()
plt.savefig("kde_bandwidth_comparison.png", dpi=120)
plt.show()

7. 演習課題（空間演習②）¶

評価対象課題です。 締め切りと提出方法は授業内で別途案内します。

課題内容¶

学内配布の「地理空間情報LLMログデータ」（llm_access_log_points.csv）を使用して、以下の成果物を作成してください。

必須成果物¶

GeoTIFF ファイル（kde_result.tif）
Silvermanの経験則で推定したバンド幅を使用すること
解像度：グリッドサイズ 300×300 以上
CRS：EPSG:4326
PNG ファイル（kde_result.png）
観測点（scatter）とKDEヒートマップを重ね合わせて表示
カラーバーを必ず付けること
タイトルに使用したバンド幅の値を明記すること

提出物一覧¶

ファイル名	説明
`kde_result.tif`	GeoTIFF形式のKDE出力
`kde_result.png`	ヒートマップ可視化（PNG）
`kde_solution.py`	作成したPythonスクリプト
`report.md`	考察（バンド幅の選択理由、密度の高い地域の解釈など）100〜200字

採点基準¶

項目	配点
GeoTIFFが正しく生成されている（CRS・解像度）	30点
PNGに観測点・ヒートマップ・カラーバーが含まれる	30点
コードが再現可能（他者の環境で動作する）	20点
考察（バンド幅の選択理由・密度パターンの解釈）	20点

ヒント¶

# バンド幅をレポートに記載するための出力例
print(f"使用したバンド幅: h = {h_silverman:.5f} 度")
print(f"（約 {h_silverman * 111:.2f} km に相当）")
# ※ 1度 ≒ 111 km（緯度方向）

まとめ¶

KDE はポイントデータから連続的な密度関数を推定する手法であり、ヒートマップ化に適している
ガウシアンカーネルが最も一般的に使われ、scikit-learnの KernelDensity や scipy の gaussian_kde で簡単に実装できる
バンド幅（h） はKDEで最重要のパラメータ。Silvermanの経験則で自動推定できるが、最終的には視覚的確認と目的に合わせた調整が必要
gdal_grid を使えばCLIからラスター生成が可能。大規模データやパイプライン処理に向いている
folium を使えばインタラクティブなWebマップとして共有できる

次回予告¶

第06回では、大規模な地理データを効率的に配信するためのベクトルタイル技術と、Tippecanoe を使った変換方法を学びます。

第05回 カーネル密度推定（KDE）¶