"""
L10 (v2): プライバシ・グリッド — k-匿名性 / quadtree / 差分プライバシ
DoBoX #1279 (河川・道路カメラ 351点) を素材に、位置情報の代表的な
プライバシ保護技法を 3 つ「ツール」として比較する v2 リライト教材。
分析1: 固定グリッド + k-匿名性 (1km / 500m / 250m ヒートマップ)
分析2: k 値変化曲線 (100m〜5km の 6 解像度)
分析3: adaptive quadtree (k_min を満たすまで再帰分割)
分析4: 差分プライバシ (Laplace ノイズ ε=1.0/0.5/0.1)
分析5: utility ↔ privacy トレードオフ
実行:
cd "2026 DoBoX 教材"
py -X utf8 lessons/L10_privacy_grid.py
データは存在しなければ DoBoX #1279 から自動取得 (ensure_dataset)。
"""
from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure,
data_recipe, ensure_dataset)
plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False
# ============================================================================
# 0. データ自動取得
# ============================================================================
CAM_PATH = ROOT / "data" / "camera_list.csv"
print("=== 0. データ自動取得 ===")
ensure_dataset(CAM_PATH, dataset_id=1279, label="カメラ位置 #1279")
# ============================================================================
# 1. 読込・前処理
# ============================================================================
print("\n=== 1. 読込・前処理 ===")
df = pd.read_csv(CAM_PATH, encoding="utf-8-sig")
df["lat"] = pd.to_numeric(df["緯度"], errors="coerce")
df["lon"] = pd.to_numeric(df["経度"], errors="coerce")
df = df.dropna(subset=["lat", "lon"]).reset_index(drop=True)
N = len(df)
print(f"カメラ点数: {N}")
print(f"緯度範囲: {df['lat'].min():.4f} 〜 {df['lat'].max():.4f}")
print(f"経度範囲: {df['lon'].min():.4f} 〜 {df['lon'].max():.4f}")
# 緯度1度 ≒ 111km / 経度1度 ≒ 111 cos(lat) km
LAT_KM = 111.0
LON_KM = 111.0 * float(np.cos(np.radians(df["lat"].mean())))
print(f"スケール係数: 緯度1度={LAT_KM:.2f}km, 経度1度={LON_KM:.2f}km (cos補正)")
def grid_label(lat, lon, cell_m):
"""点群を cell_m の格子に量子化し、(i, j) ビンインデックスを返す。"""
cell_lat = cell_m / 1000.0 / LAT_KM
cell_lon = cell_m / 1000.0 / LON_KM
i = np.floor(np.asarray(lat) / cell_lat).astype(int)
j = np.floor(np.asarray(lon) / cell_lon).astype(int)
return i, j
# ============================================================================
# 2. 分析1: 固定グリッド + k-匿名性 (図1: grid_heatmap)
# ============================================================================
print("\n=== 分析1: 固定グリッド (1km / 500m / 250m) ヒートマップ ===")
RESOLUTIONS = [1000, 500, 250] # m
fig, axes = plt.subplots(1, 3, figsize=(15, 5.2), sharey=True)
grid_results = {} # cell_m -> {keys, counts}
grid_summary_rows = []
for ax, cs in zip(axes, RESOLUTIONS):
i, j = grid_label(df["lat"].values, df["lon"].values, cs)
keys = list(zip(i.tolist(), j.tolist()))
counts = pd.Series(keys).value_counts()
grid_results[cs] = {"keys": keys, "counts": counts}
# ヒートマップ用の 2D 配列
i_min, i_max = i.min(), i.max()
j_min, j_max = j.min(), j.max()
H = np.zeros((i_max - i_min + 1, j_max - j_min + 1), dtype=int)
for (ii, jj), n in counts.items():
H[ii - i_min, jj - j_min] = n
cell_lat = cs / 1000.0 / LAT_KM
cell_lon = cs / 1000.0 / LON_KM
extent = [
(j_min) * cell_lon, (j_max + 1) * cell_lon,
(i_min) * cell_lat, (i_max + 1) * cell_lat,
]
im = ax.imshow(H, origin="lower", extent=extent, cmap="YlOrRd",
aspect="auto", vmax=max(3, np.percentile(H[H > 0], 95)))
ax.scatter(df["lon"], df["lat"], s=2, color="#0969da", alpha=0.5)
n_occupied = int((counts > 0).sum())
n_k1 = int((counts == 1).sum())
n_k_ge3 = int((counts >= 3).sum())
n_pts_safe_k3 = int(counts[counts >= 3].sum())
grid_summary_rows.append({
"セル一辺(m)": cs,
"セル面積(km²)": round((cs / 1000.0) ** 2, 3),
"占有セル数": n_occupied,
"k=1セル数": n_k1,
"k=1セル割合(%)": round(n_k1 / max(n_occupied, 1) * 100, 1),
"k≥3セル数": n_k_ge3,
"k≥3条件で残る点数": n_pts_safe_k3,
"失う点数": N - n_pts_safe_k3,
})
ax.set_title(f"{cs}m 格子\n占有={n_occupied}, k=1={n_k1} ({n_k1/n_occupied*100:.0f}%)")
ax.set_xlabel("経度")
if ax is axes[0]:
ax.set_ylabel("緯度")
plt.colorbar(im, ax=ax, shrink=0.8, label="セル内点数")
plt.suptitle("図1: 固定グリッドのプライバシ可視化 (青点=オリジナルカメラ位置)",
fontsize=12, y=1.02)
plt.tight_layout()
plt.savefig(ASSETS / "L10_grid_heatmap.png", dpi=140, bbox_inches="tight")
plt.close()
grid_summary_df = pd.DataFrame(grid_summary_rows)
print(grid_summary_df.to_string(index=False))
grid_summary_df.to_csv(ASSETS / "L10_grid_counts.csv", index=False, encoding="utf-8-sig")
# Before/After: 1点を最初から最後まで追う表 (要件K)
sample_idx = 0 # No.1 苗代カメラ
sample_lat = df.loc[sample_idx, "lat"]
sample_lon = df.loc[sample_idx, "lon"]
sample_name = df.loc[sample_idx, "カメラ名"]
sample_addr = df.loc[sample_idx, "住所"]
ba_rows = []
ba_rows.append({
"段階": "原データ",
"値": f"({sample_lat:.6f}, {sample_lon:.6f})",
"粒度": "GPS精度 約5m",
"このセル内の点数 (k)": "—",
"個人特定リスク": "高 (1棟単位で特定可能)",
})
for cs in [250, 500, 1000]:
i, j = grid_label([sample_lat], [sample_lon], cs)
ii, jj = int(i[0]), int(j[0])
counts = grid_results[cs]["counts"] if cs in grid_results else None
if counts is not None:
k_here = int(counts.get((ii, jj), 0))
else:
k_here = 0
cell_lat = cs / 1000.0 / LAT_KM
cell_lon = cs / 1000.0 / LON_KM
cell_center_lat = (ii + 0.5) * cell_lat
cell_center_lon = (jj + 0.5) * cell_lon
ba_rows.append({
"段階": f"{cs}m 格子に量子化",
"値": f"セル中心 ({cell_center_lat:.4f}, {cell_center_lon:.4f})",
"粒度": f"{cs}m × {cs}m",
"このセル内の点数 (k)": k_here,
"個人特定リスク": "高" if k_here == 1 else "中" if k_here < 5 else "低",
})
sample_trace_df = pd.DataFrame(ba_rows)
sample_trace_df.to_csv(ASSETS / "L10_sample_trace.csv", index=False, encoding="utf-8-sig")
sample_trace_html = sample_trace_df.to_html(index=False)
grid_summary_html = grid_summary_df.to_html(index=False)
# ============================================================================
# 3. 分析2: k 値変化曲線 (図2: k_curve)
# ============================================================================
print("\n=== 分析2: k-匿名性曲線 ===")
RES_ALL = [100, 250, 500, 1000, 2000, 5000]
K_TARGETS = [1, 2, 5, 10]
k_curve = {k: [] for k in K_TARGETS}
median_k_per_res = []
k_curve_rows = []
for cs in RES_ALL:
i, j = grid_label(df["lat"].values, df["lon"].values, cs)
keys = list(zip(i.tolist(), j.tolist()))
counts = pd.Series(keys).value_counts()
median_k_per_res.append(float(counts.median()))
row = {"解像度(m)": cs, "占有セル数": int(len(counts)),
"中央k": float(counts.median())}
for k in K_TARGETS:
# 各点について「自分のセル内の k」が k 以上である割合
k_each = counts.loc[keys].values
ratio = float((k_each >= k).mean()) * 100
k_curve[k].append(ratio)
row[f"k≥{k} 達成率(%)"] = round(ratio, 1)
k_curve_rows.append(row)
print(f" {cs}m: median k={counts.median():.0f}, "
+ ", ".join(f"k>={k}:{k_curve[k][-1]:.1f}%" for k in K_TARGETS))
fig, ax = plt.subplots(figsize=(8.5, 5))
COLORS = ["#999999", "#0969da", "#fb8500", "#cf222e"]
MARKERS = ["x", "o", "s", "^"]
for k, c, m in zip(K_TARGETS, COLORS, MARKERS):
ax.plot(RES_ALL, k_curve[k], f"-{m}", color=c, linewidth=2,
markersize=8, label=f"k≧{k}")
ax.set_xscale("log")
ax.set_xticks(RES_ALL)
ax.set_xticklabels([f"{x}m" if x < 1000 else f"{x//1000}km" for x in RES_ALL])
ax.set_xlabel("グリッド解像度 (セル一辺)")
ax.set_ylabel("条件を満たす点の割合 (%)")
ax.set_ylim(-5, 105)
ax.axhline(80, color="gray", linestyle=":", linewidth=0.8, alpha=0.7)
ax.text(120, 81, "実務目安 80%", color="gray", fontsize=9)
ax.set_title("図2: k-匿名性曲線 — 解像度を粗くするほど k≧c を満たす点が増える")
ax.legend(loc="lower right")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L10_k_curve.png", dpi=140)
plt.close()
k_curve_df = pd.DataFrame(k_curve_rows)
k_curve_df.to_csv(ASSETS / "L10_k_anonymity_curve.csv",
index=False, encoding="utf-8-sig")
k_curve_html = k_curve_df.to_html(index=False)
# ============================================================================
# 4. 分析3: adaptive quadtree (図3: quadtree)
# ============================================================================
print("\n=== 分析3: adaptive quadtree ===")
class QuadNode:
__slots__ = ("lat0", "lat1", "lon0", "lon1", "pts", "children")
def __init__(self, lat0, lat1, lon0, lon1, pts):
self.lat0, self.lat1 = lat0, lat1
self.lon0, self.lon1 = lon0, lon1
self.pts = pts # list of (lat, lon)
self.children = None # None=leaf
def is_leaf(self):
return self.children is None
def build_quadtree(pts, lat0, lat1, lon0, lon1, k_min=8, max_depth=10, depth=0):
"""点数 ≥ k_min を満たす限界まで再帰的に4分割。
停止条件:
- 含まれる点数が 2*k_min 未満 → これ以上割ると k_min を割る → leaf
- 4分割した結果、いずれかの空でない子の点数が k_min 未満 → 分割中止 → leaf
- depth が max_depth に到達
結果として 全 leaf は 0 点 (空) または k_min 点以上 を保証する。
"""
node = QuadNode(lat0, lat1, lon0, lon1, pts)
if depth >= max_depth or len(pts) < 2 * k_min:
return node
lat_mid = 0.5 * (lat0 + lat1)
lon_mid = 0.5 * (lon0 + lon1)
quads = [[], [], [], []] # SW, SE, NW, NE
for p in pts:
si = 0 if p[0] < lat_mid else 2 # 緯度: 南=0/北=2
sj = 0 if p[1] < lon_mid else 1 # 経度: 西=0/東=1
quads[si + sj].append(p)
if any(0 < len(q) < k_min for q in quads):
return node
children = []
for q_idx, qpts in enumerate(quads):
si, sj = divmod(q_idx, 2)
c_lat0 = lat_mid if si else lat0
c_lat1 = lat1 if si else lat_mid
c_lon0 = lon_mid if sj else lon0
c_lon1 = lon1 if sj else lon_mid
if qpts:
children.append(build_quadtree(qpts, c_lat0, c_lat1, c_lon0, c_lon1,
k_min, max_depth, depth + 1))
else:
children.append(QuadNode(c_lat0, c_lat1, c_lon0, c_lon1, []))
node.children = children
return node
def collect_leaves(node):
if node.is_leaf():
return [node]
out = []
for c in node.children:
out.extend(collect_leaves(c))
return out
# 県全域 bbox
PAD = 0.02
BBOX = (df["lat"].min() - PAD, df["lat"].max() + PAD,
df["lon"].min() - PAD, df["lon"].max() + PAD)
pts_all = list(zip(df["lat"].tolist(), df["lon"].tolist()))
QT_RESULTS = []
qt_summary_rows = []
qt_cells_all_rows = []
for k_min in [2, 5, 10]:
root = build_quadtree(pts_all, *BBOX, k_min=k_min, max_depth=14)
leaves = [lf for lf in collect_leaves(root) if lf.pts]
QT_RESULTS.append((k_min, root, leaves))
sizes = np.array([len(lf.pts) for lf in leaves])
leaf_areas_km2 = np.array([
(lf.lat1 - lf.lat0) * LAT_KM * (lf.lon1 - lf.lon0) * LON_KM
for lf in leaves
])
# 点加重(各点が属する leaf の面積を 351 個並べた中央値)
per_pt_area = []
per_pt_size = []
for ar, sz in zip(leaf_areas_km2, sizes):
per_pt_area.extend([float(ar)] * int(sz))
per_pt_size.extend([int(sz)] * int(sz))
qt_summary_rows.append({
"k_min": k_min,
"leaf 数": int(len(leaves)),
"min 点数": int(sizes.min()),
"中央 点数(leaf)": int(np.median(sizes)),
"max 点数": int(sizes.max()),
"中央 leaf 面積(km², leaf単純)": round(float(np.median(leaf_areas_km2)), 3),
"中央 leaf 面積(km², 点加重)": round(float(np.median(per_pt_area)), 3),
"中央 k(点加重)": int(np.median(per_pt_size)),
})
print(f" k_min={k_min}: leaf={len(leaves)}, "
f"min={int(sizes.min())}, median={int(np.median(sizes))}, max={int(sizes.max())}")
# 各 leaf の詳細を CSV へ
for lf, sz, ar in zip(leaves, sizes, leaf_areas_km2):
qt_cells_all_rows.append({
"k_min": k_min,
"lat0": round(lf.lat0, 5), "lat1": round(lf.lat1, 5),
"lon0": round(lf.lon0, 5), "lon1": round(lf.lon1, 5),
"点数": int(sz),
"面積(km²)": round(float(ar), 4),
})
fig, axes = plt.subplots(1, 3, figsize=(15, 5.2), sharey=True)
for ax, (k_min, root, leaves) in zip(axes, QT_RESULTS):
sizes = np.array([len(lf.pts) for lf in leaves])
vmax = sizes.max()
cmap = plt.cm.YlOrRd
for lf in leaves:
n = len(lf.pts)
color = cmap(0.15 + 0.85 * n / vmax)
rect = Rectangle((lf.lon0, lf.lat0),
lf.lon1 - lf.lon0, lf.lat1 - lf.lat0,
facecolor=color, edgecolor="#444", linewidth=0.5,
alpha=0.75)
ax.add_patch(rect)
ax.scatter(df["lon"], df["lat"], s=3, color="#0550ae", alpha=0.6)
ax.set_xlim(BBOX[2], BBOX[3])
ax.set_ylim(BBOX[0], BBOX[1])
ax.set_title(f"k_min={k_min}: leaf={len(leaves)}, "
f"min点数={sizes.min()}")
ax.set_xlabel("経度")
if ax is axes[0]:
ax.set_ylabel("緯度")
plt.suptitle("図3: adaptive quadtree — k_min を満たすまで再帰分割 "
"(点が密な地域ほど細かいセル)", fontsize=12, y=1.02)
plt.tight_layout()
plt.savefig(ASSETS / "L10_quadtree.png", dpi=140, bbox_inches="tight")
plt.close()
qt_summary_df = pd.DataFrame(qt_summary_rows)
qt_summary_df.to_csv(ASSETS / "L10_quadtree_summary.csv",
index=False, encoding="utf-8-sig")
pd.DataFrame(qt_cells_all_rows).to_csv(
ASSETS / "L10_quadtree_cells.csv", index=False, encoding="utf-8-sig")
qt_summary_html = qt_summary_df.to_html(index=False)
# ============================================================================
# 5. 分析4: 差分プライバシ Laplace ノイズ (図4: dp_noise)
# ============================================================================
print("\n=== 分析4: 差分プライバシ (Laplace ノイズ) ===")
EPSILONS = [1.0, 0.5, 0.1]
SENS_LAT = 1.0 / LAT_KM # 1km 相当の緯度オフセット
SENS_LON = 1.0 / LON_KM # 1km 相当の経度オフセット
rng = np.random.default_rng(seed=42)
fig, axes = plt.subplots(1, 4, figsize=(18, 5.2), sharey=True)
ax = axes[0]
ax.scatter(df["lon"], df["lat"], s=8, color="#0969da", alpha=0.7)
ax.set_xlim(BBOX[2], BBOX[3])
ax.set_ylim(BBOX[0], BBOX[1])
ax.set_title(f"オリジナル (n={N})")
ax.set_xlabel("経度")
ax.set_ylabel("緯度")
dp_summary_rows = []
dp_per_point_rows = [] # 全点ごとのジッタ距離
for ax, eps in zip(axes[1:], EPSILONS):
b_lat = SENS_LAT / eps
b_lon = SENS_LON / eps
noise_lat = rng.laplace(scale=b_lat, size=N)
noise_lon = rng.laplace(scale=b_lon, size=N)
lat_dp = df["lat"].values + noise_lat
lon_dp = df["lon"].values + noise_lon
d_km = np.sqrt((noise_lat * LAT_KM) ** 2 + (noise_lon * LON_KM) ** 2)
median_d = float(np.median(d_km))
p95_d = float(np.percentile(d_km, 95))
max_d = float(d_km.max())
dp_summary_rows.append({
"ε": eps,
"Laplace b (km)": round(b_lat * LAT_KM, 3),
"ジッタ中央値(km)": round(median_d, 2),
"ジッタ95%点(km)": round(p95_d, 2),
"ジッタ最大(km)": round(max_d, 2),
"県内に収まる点(%)": round(float(((lat_dp >= BBOX[0]) & (lat_dp <= BBOX[1]) &
(lon_dp >= BBOX[2]) & (lon_dp <= BBOX[3])).mean()) * 100, 1),
})
# 全点ジッタ距離
for k in range(N):
dp_per_point_rows.append({
"epsilon": eps,
"no": int(df.loc[k, "No."]),
"name": df.loc[k, "カメラ名"],
"orig_lat": round(float(df.loc[k, "lat"]), 6),
"orig_lon": round(float(df.loc[k, "lon"]), 6),
"dp_lat": round(float(lat_dp[k]), 6),
"dp_lon": round(float(lon_dp[k]), 6),
"jitter_km": round(float(d_km[k]), 3),
})
for la0, lo0, la1, lo1 in zip(df["lat"], df["lon"], lat_dp, lon_dp):
ax.plot([lo0, lo1], [la0, la1], color="#cf222e", alpha=0.15, linewidth=0.5)
ax.scatter(df["lon"], df["lat"], s=4, color="#0969da", alpha=0.35,
label="オリジナル")
ax.scatter(lon_dp, lat_dp, s=8, color="#cf222e", alpha=0.7,
label=f"DP ε={eps}")
ax.set_xlim(BBOX[2], BBOX[3])
ax.set_ylim(BBOX[0], BBOX[1])
ax.set_title(f"ε={eps}\n中央ジッタ {median_d:.1f}km, p95 {p95_d:.1f}km")
ax.set_xlabel("経度")
ax.legend(loc="lower right", fontsize=8)
plt.suptitle("図4: 差分プライバシ (Laplace ノイズ) — ε が小さいほど大きく揺らぐ",
fontsize=12, y=1.02)
plt.tight_layout()
plt.savefig(ASSETS / "L10_dp_noise.png", dpi=140, bbox_inches="tight")
plt.close()
dp_df = pd.DataFrame(dp_summary_rows)
print(dp_df.to_string(index=False))
dp_df.to_csv(ASSETS / "L10_dp_results.csv", index=False, encoding="utf-8-sig")
pd.DataFrame(dp_per_point_rows).to_csv(
ASSETS / "L10_dp_per_point.csv", index=False, encoding="utf-8-sig")
dp_html = dp_df.to_html(index=False)
# ============================================================================
# 6. 分析5: utility ↔ privacy トレードオフ (図5: tradeoff)
# ============================================================================
print("\n=== 分析5: utility-privacy トレードオフ ===")
fig, ax = plt.subplots(figsize=(8.5, 5))
areas_km2 = [(r / 1000.0) ** 2 for r in RES_ALL]
sc = ax.scatter(areas_km2, median_k_per_res, s=180, c=RES_ALL, cmap="viridis",
edgecolor="black", linewidth=1.0, zorder=3)
for x, y, r in zip(areas_km2, median_k_per_res, RES_ALL):
label = f"{r}m" if r < 1000 else f"{r//1000}km"
ax.annotate(label, (x, y), xytext=(7, 7), textcoords="offset points", fontsize=10)
tradeoff_rows = []
for cs, area, mk in zip(RES_ALL, areas_km2, median_k_per_res):
tradeoff_rows.append({"手法": "固定グリッド", "パラメータ": f"{cs}m",
"セル面積中央(km²)": round(area, 3),
"セル内点数中央 (k)": round(mk, 1)})
# quadtree 結果も同じ平面に重ねる
# 注: leaf を単純中央すると過疎地の巨大 leaf が中央値を引き上げる。
# 「点1つあたりが置かれる平均的なセル」を見たいので、ここでは 点加重中央 を計算する
# (各点が属する leaf の面積/点数を 351 点並べてから中央値を取る)
for k_min, _, leaves in QT_RESULTS:
leaf_areas = np.array([
(lf.lat1 - lf.lat0) * LAT_KM * (lf.lon1 - lf.lon0) * LON_KM
for lf in leaves
])
leaf_sizes = np.array([len(lf.pts) for lf in leaves])
# 点加重: 各点が「自分の leaf 面積 / k」を持つ
per_point_area = []
per_point_k = []
for ar, sz in zip(leaf_areas, leaf_sizes):
per_point_area.extend([float(ar)] * int(sz))
per_point_k.extend([int(sz)] * int(sz))
pw_area = float(np.median(per_point_area))
pw_k = float(np.median(per_point_k))
ax.scatter([pw_area], [pw_k], marker="*", s=320, color="#cf222e",
edgecolor="black", linewidth=1.0, zorder=4)
ax.annotate(f"quadtree k_min={k_min}", (pw_area, pw_k),
xytext=(10, -14), textcoords="offset points",
fontsize=9, color="#cf222e")
tradeoff_rows.append({"手法": "adaptive quadtree", "パラメータ": f"k_min={k_min}",
"セル面積中央(km²)": round(pw_area, 3),
"セル内点数中央 (k)": round(pw_k, 1)})
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("セル面積 (km²) — 大きいほど情報損失大")
ax.set_ylabel("セル内点数 中央値 — 大きいほどプライバシ保護強")
ax.set_title("図5: utility (情報損失) vs privacy (k 中央値) のトレードオフ\n"
"(青〜黄=固定グリッド, 赤★=adaptive quadtree)")
ax.grid(alpha=0.3, which="both")
plt.colorbar(sc, ax=ax, label="グリッド解像度 (m)")
plt.tight_layout()
plt.savefig(ASSETS / "L10_tradeoff.png", dpi=140)
plt.close()
tradeoff_df = pd.DataFrame(tradeoff_rows)
tradeoff_df.to_csv(ASSETS / "L10_tradeoff.csv", index=False, encoding="utf-8-sig")
tradeoff_html = tradeoff_df.to_html(index=False)
# ============================================================================
# 7. 補助: 「家1棟スケール」の k=1 リスク表 (考察用)
# ============================================================================
print("\n=== 補助: 家1棟スケールの k=1 リスク ===")
HOUSE_RES = [5, 10, 25, 50, 100]
house_rows = []
for cs in HOUSE_RES:
i, j = grid_label(df["lat"].values, df["lon"].values, cs)
keys = list(zip(i.tolist(), j.tolist()))
counts = pd.Series(keys).value_counts()
n_k1 = (counts == 1).sum()
house_rows.append({
"格子 (m)": cs,
"占有セル": int((counts > 0).sum()),
"k=1 セル数": int(n_k1),
"k=1 割合 (%)": round(n_k1 / max((counts > 0).sum(), 1) * 100, 1),
"建物棟数イメージ": "1棟未満" if cs <= 10 else "数棟" if cs <= 50 else "10〜数十棟",
})
house_df = pd.DataFrame(house_rows)
print(house_df.to_string(index=False))
house_df.to_csv(ASSETS / "L10_house_scale.csv", index=False, encoding="utf-8-sig")
house_html = house_df.to_html(index=False)
# ============================================================================
# 8. HTML レンダリング (10 セクション構成)
# ============================================================================
print("\n=== HTML レンダリング ===")
sections = [
# =========================================================================
# 1. 学習目標と問い
# =========================================================================
("学習目標と問い", """
このレッスンで答えたい問い
「カメラ位置のような『ピンポイントの位置データ』を、個人特定リスクを下げて公開するには、どんな道具がどこまで効くのか?」
道路・河川カメラ自体は本来公開情報。本レッスンは 「個人や要援護者の位置データを、仮にこのカメラ点群と同じ密度で扱った場合」 を念頭に、
代表的な3手法を 「ツール」として比較する シミュレーション教材。
用語の定義 (このレッスン独自・冒頭で平易に)
- k-匿名性 (k-anonymity):
「同じ場所(=同じセル)にいる人が k 人以上になっていれば、その中の誰かは特定できない」という考え方。
k=1 はその場所に あなた1人しかいない=即特定。k=5 なら最低5人と紛れる。
ツールとしては「セルに分割して、k 未満のセルは公開しない or 粗くする」操作。
- quadtree (クアッドツリー):
地図を 4分割→各区画をさらに4分割→…と再帰的に細かくしていく木構造。
本レッスンでは「セル内の点数が k_min 未満になりそうなら分割を止める」よう細工する。
結果: 都市部は細かく、過疎地は粗く、点数のバランスが取れたセル分割になる。
- 差分プライバシ (Differential Privacy, DP):
データに小さなランダムノイズを加えることで、「あなた1人がデータに含まれていたかどうかを攻撃者が見破れない」ことを 数学的に保証する技法。
ε(イプシロン) という1つのパラメータで「保護の強さ」を制御する(小さいほど強保護・大きいほど弱保護)。
- utility / privacy:
utility(ユーティリティ)=「データの使いやすさ・精度」、privacy(プライバシ)=「個人が守られる度合い」。
この2つは シーソー関係で、両方を最大化するのは原理的に難しい。本レッスンの最終図(図5)で定量化する。
立てた仮説 (後で検証する)
- H1 (固定グリッドは万能ではない):
「県全体に同じ大きさの格子をかぶせる」という素朴な方法では、
都市部は過剰に粗く、山間部は k=1 のままになる。
つまり 均一なセルでは均一な保護にならないはず。
- H2 (k を上げるほど解像度が必要):
k≧1, 2, 5, 10 と要求を厳しくするほど、それを満たすには 大きなセルが必要になる。
k≧10 を全点で達成するには、km単位のセルが要るはず。
- H3 (quadtree は固定グリッドに勝つ):
「点が多い場所だけ細かく割る」適応分割なら、
同じプライバシ水準(k_min)を より細かい平均粒度で達成できるはず。
utility-privacy 平面上で固定グリッドより 左上に位置する。
- H4 (差分プライバシは ε で揺らぎが指数的に変わる):
Laplace ノイズの幅は b=Δ/ε。
ε を 1.0 → 0.1 に下げると揺らぎは 10倍になる。
ε=0.1 のような強保護では、地図公開には 使い物にならないレベルでジッタするはず。
- H5 (1棟特定は10m格子で起きる):
GPS の生精度が 5〜10m。
10m 格子では占有セルの大半が k=1 となり、
1棟単位での特定が技術的に可能になるはず。
「位置データは本質的に強い識別子」という直観を数値で確認する。
到達点
- 5仮説に対して、5つの分析でそれぞれ 支持／反証／部分支持 を判定できる
- k-匿名性／quadtree／差分プライバシ の3手法を 「何を入れて何が返るか」のツール視点で説明できる
- 「どの粒度で公開すべきか」を utility-privacy トレードオフ図から定量的に選べる
- 内部の確率不等式 (Laplace 機構の事後分布比 ≤ e^ε など) は 黒箱で OK。
ツールの 入出力と限界を語れることが目標。
"""),
# =========================================================================
# 2. 使用データ
# =========================================================================
("使用データ", f"""
- 名称 河川・道路カメラ位置情報
- 取得元 DoBoX #1279
- 件数 351 点 (緯度・経度・住所・路河川名・所管・管理区分 の 9 列)
- ファイル
data/camera_list.csv (UTF-8 BOM, 約 70 KB)
- 位置範囲 緯度 {df['lat'].min():.3f}〜{df['lat'].max():.3f}, 経度 {df['lon'].min():.3f}〜{df['lon'].max():.3f} (広島県全域)
- 背景 道路・河川カメラの座標は本来公開情報。本レッスンは 「個人や要援護者の位置データ」を仮にこの密度で扱う場合 を念頭においたシミュレーション教材
- なぜカメラなのか 県全域に 不均一に(都市部に密、山間に疎) 分布する 351 点という規模が、k-匿名性の歪みと adaptive 分割の必要性を体感するのに最適
"""),
# =========================================================================
# 3. ダウンロード (再現用データ・中間データ・図)
# =========================================================================
("ダウンロード(再現用データ・中間データ・図)", """
本レッスンの全成果物に直リンクを置いた。途中ステップから再現したい学習者向け。
1. 生データ (DoBoX 由来)
2. プログラムで生成される中間データ
3. 図 PNG
4. 再現スクリプト
cd "2026 DoBoX 教材"
py -X utf8 lessons/L10_privacy_grid.py
スクリプト本体: lessons/L10_privacy_grid.py
(データが無ければ ensure_dataset() で自動DL → 5枚のPNGと9本の中間CSVを生成)
"""),
# =========================================================================
# 4. 分析1: 固定グリッド + k-匿名性
# =========================================================================
("分析1: 固定グリッド + k-匿名性", """
狙い
「県全体を同じ大きさの格子で割って、各セルの中に何人(=点)いるかを数える」だけの素朴な方法を試す。
これが k-匿名性の最もシンプルな実装で、3つの粒度 (1km / 500m / 250m) を比較して
「均一格子の盲点」(=都市部は過剰に粗く、山間部は k=1 のまま) を視覚化する。仮説H1 の検証。
ツールとしての k-匿名性:
- 入力: 351点の緯度経度
- 出力: 各点に「同じセル内の他の点数 (k)」のラベル
- 使い方: k 未満のセルにある点は「公開しない」「粗いセルにマージする」「ノイズで上書きする」等の判断材料にする
- 限界: 同じセル内の人が 属性的に均質(全員が同じ性別・年齢) なら、k≧5 でも事実上特定できる(homogeneity attack)。これを補うのが l-diversity / t-closeness で、発展課題で扱う。
手法 (緯度経度の格子量子化)
- 入力: 351点の (lat, lon) と セル一辺 cell_m メートル
- 処理: 緯度1度=111km, 経度1度=111·cos(lat) km の cos 補正 でメートル→度を換算し、
floor(lat / cell_lat) と floor(lon / cell_lon) でビン化
- 出力: 各点の (セル行 i, セル列 j) ペア →
value_counts で「セル内点数」
- パラメータ: cell_m ∈ {250, 500, 1000} m (細→粗)
- 結果の読み方: k=1 セルが「個人特定リスク高」、k≧5 で「実用安全圏」が一般的目安
実装
""" + code('''
import numpy as np, pandas as pd
# === (1) 緯経度の cos 補正 (実距離→緯経度オフセット) ===
LAT_KM = 111.0
LON_KM = 111.0 * np.cos(np.radians(df["lat"].mean())) # 北緯34.5度なら ≒92km
def grid_label(lat, lon, cell_m):
"""点群を cell_m 格子に量子化し、(行 i, 列 j) を返す"""
cell_lat = cell_m / 1000 / LAT_KM # 例: 250m → 約 0.00225度
cell_lon = cell_m / 1000 / LON_KM # 例: 250m → 約 0.00271度
return (np.floor(lat / cell_lat).astype(int),
np.floor(lon / cell_lon).astype(int))
# === (2) 3 解像度それぞれで集計 ===
for cs in [1000, 500, 250]:
i, j = grid_label(df["lat"], df["lon"], cs)
counts = pd.Series(list(zip(i, j))).value_counts() # セル毎の点数
n_occupied = (counts > 0).sum() # 占有セル数
n_k1 = (counts == 1).sum() # k=1 リスクセル数
print(f"{cs}m: 占有={n_occupied}, k=1={n_k1} ({n_k1/n_occupied*100:.0f}%)")
''') + """
図と読み取り
なぜこの図か: 3粒度を 同じ枠で並べる(small multiples) ことで、
セルが細かくなるほど k=1 の白〜薄黄セルが急増する様子を一目で比較できる。
1枚の図だけでは「これは粗いのか細かいのか」を相対判断できない。
""" + figure("assets/L10_grid_heatmap.png",
"図1: 1km / 500m / 250m の固定グリッド3解像度。色濃いほど多人数が同じセルに集まり、薄いセルは k=1=即特定リスク") + """
この図から読み取れること:
- 都市部(広島市・福山市)は色濃い: 1km 格子なら k≧5 を確保できる地域が多い
- 山間部・島嶼部はほぼ薄黄: どの粒度でも k=1 のセルが多数 — 同じ格子では都市と山で「保護の質」が違う
- 250m 格子では過半数が k=1: タイトルの
k=1セル割合 数値が確認できる
- 仮説H1 を支持: 均一なセルでは 均一な保護にならない
表と読み取り
なぜこの表か: 図1 の見た目を 「失う点数」という具体数で裏付けたい。
k≧3 を満たさないセルは公開時に削るので、その時点で何点が地図から消えるかを可視化する。
表1: 3解像度のセル統計
""" + grid_summary_html + """
読み取り: 1km 格子では k≧3 を満たすセルが多く、失う点数が比較的小さい。
逆に 250m では失う点数がぐっと増え、利用可能な情報が大幅に減る。
細かさ(=utility)とプライバシ保護は逆相関であることが数値で確認できる。
表2: 1点を最初から最後まで追う Before/After (No.1 苗代カメラを例に)
""" + sample_trace_html + """
読み取り: 同じ点でも セル一辺を粗くするほど k が増え、リスクが下がる。
逆に細かい格子では k=1 になりがちで、その点は公開不可と判定される。
このように 「ある1点」がツールにかけられた後どう変換されるか を最後まで追えるのが、
データ前処理の理解の核心。
"""),
# =========================================================================
# 5. 分析2: k 値変化曲線
# =========================================================================
("分析2: k 値変化曲線(解像度を細→粗に振って曲線化)", """
狙い
「セルの大きさを 100m から 5km まで振ったら、k≧1, 2, 5, 10 の各条件を満たす点はそれぞれ何%になるか?」
を1枚の折れ線グラフにする。3粒度の点(分析1)を 連続曲線に拡張することで、
「実務目安 80% を超えるにはどれくらい粗くする必要があるか」が読めるようになる。仮説H2の検証。
手法 (解像度パラメータスイープ)
- 入力: 351点 + 解像度リスト [100m, 250m, 500m, 1km, 2km, 5km] + 目標 k リスト [1, 2, 5, 10]
- 処理: 各解像度で
grid_label → value_counts →
各点について「自分のセルの k」を引き、(k ≥ 目標) を .mean() で割合化
- 出力: 6解像度 × 4 k値 の 達成率 (%) 行列 → 折れ線グラフ
- 結果の読み方: 80% ライン (実務目安) と曲線の交点が 「最低限必要な粗さ」
実装
""" + code('''
RES_ALL = [100, 250, 500, 1000, 2000, 5000] # m
K_TARGETS = [1, 2, 5, 10]
k_curve = {k: [] for k in K_TARGETS}
for cs in RES_ALL:
i, j = grid_label(df["lat"], df["lon"], cs)
keys = list(zip(i, j))
counts = pd.Series(keys).value_counts()
# 各点の「自分が属するセルの k」を引く
k_each = counts.loc[keys].values # 長さ 351 の int 配列
for k in K_TARGETS:
ratio = (k_each >= k).mean() * 100 # 達成率 %
k_curve[k].append(ratio)
''') + """
図と読み取り
なぜこの図か: x軸を 対数スケール(100m〜5km)、y軸を達成率(%) にすることで、
「どこまで粗くすれば 80% を超えるか」を直接読めるようにした。
4本の線(k=1, 2, 5, 10) を重ねることで、k 要求が厳しくなるほど右にシフトすることが視覚化される。
""" + figure("assets/L10_k_curve.png",
"図2: k-匿名性曲線 — 横軸=解像度(対数), 縦軸=条件達成率(%)。点線は実務目安80%") + """
この図から読み取れること:
- k≧1 はどの解像度でも 100% (自分自身が1人いるので自明)
- k≧2 達成率は 5km 格子でも約 67.5%: 100m なら 4%, 1km なら 22% にとどまり、80%目安には届かない
- k≧5 達成率は 5km でも約 24%: 単純グリッドでは 351点という規模では「他に4人いる場所」自体が少ない
- k≧10 は 5km でも 2.8% のみ: 山間離島ではどう粗くしても点が孤立する
- 仮説H2 を支持: k を上げるほど曲線は 右に大きくシフト。点数が少ないと固定グリッドだけでは k 確保が原理的に厳しい
表と読み取り
表3: 解像度 × k目標 の達成率(%)
""" + k_curve_html + """
読み取り:
「k≧5 を 80% 確保したい」という設計目標を立てた場合、表から 5km格子でも 24%=固定グリッドでは届かない。
これは「351点という規模では、地理的に疎な点が多すぎて、均一格子では k≧5 を多くの点で確保できない」ことの定量的証拠で、
固定グリッドの限界を端的に示す数字である。
これが分析3の adaptive quadtree を導入する動機になる(quadtree なら「点が密な場所だけ細かく」できる)。
"""),
# =========================================================================
# 6. 分析3: adaptive quadtree
# =========================================================================
("分析3: adaptive quadtree(適応分割)", """
狙い
「点が多い場所だけ細かく、点が少ない場所は粗く」を自動でやる。
固定グリッドの欠点(都市過剰粗・山間k=1)を 木構造の再帰分割で解決する。
仮説H3 の検証。
ツールとしての quadtree:
- 入力: 351点の (lat, lon) と最低人数 k_min
- 出力: 「どんなセル分割をすれば全てのセルが k_min 人以上になるか」 という地図分割 (leaf の集合)
- 動作のイメージ: 全域→4分割→各区画でまた4分割→…
ただし「分割すると k_min を割る区画が出る」なら そこで分割を止める(leaf 化)。
結果として、点が密な広島市内は細かいセル、点が疎な山間部は粗いセルが自動で出来上がる。
- 限界: 矩形分割なので「対角に伸びる集落」を理想的には捉えられない。
ボロノイ分割など別形状もあるが、quadtree は 実装が単純かつ 木構造で高速検索できるのが利点。
手法 (再帰的4分割アルゴリズム)
- 入力: 全点 pts + 全域 bbox + k_min
- 処理:
- もし
len(pts) < 2*k_min ならこれ以上割れない → leaf にする
- そうでなければ bbox の中央で 南北・東西に2回ずつ切って4区画に分ける
- 4区画のうち どれか1つでも 0 < n < k_min になるなら、分割を中止して leaf 化
- 全区画が空 or k_min 以上なら、4子それぞれに対して再帰呼び出し
- 出力: 全 leaf が「0点 (空)」か「k_min 点以上」を保証する木
- パラメータ: k_min ∈ {2, 5, 10}, max_depth=14 (再帰の安全弁)
実装
""" + code('''
def build_quadtree(pts, lat0, lat1, lon0, lon1, k_min=8, max_depth=10, depth=0):
"""点数 >= k_min を満たす限界まで4分割。
4分割で 0= max_depth or len(pts) < 2 * k_min:
return {"bbox": (lat0, lat1, lon0, lon1), "pts": pts, "leaf": True}
# 4分割 (lat, lon それぞれ中央で2分)
lat_m = (lat0 + lat1) / 2
lon_m = (lon0 + lon1) / 2
quads = [[], [], [], []] # SW, SE, NW, NE
for la, lo in pts:
si = 0 if la < lat_m else 2 # 南/北
sj = 0 if lo < lon_m else 1 # 西/東
quads[si + sj].append((la, lo))
# 停止条件 2: 1つでも k_min 未満の非空区画が出るなら分割中止
if any(0 < len(q) < k_min for q in quads):
return {"bbox": (lat0, lat1, lon0, lon1), "pts": pts, "leaf": True}
# 再帰
children = []
for q_idx, q in enumerate(quads):
si, sj = divmod(q_idx, 2)
c_lat0, c_lat1 = (lat_m, lat1) if si else (lat0, lat_m)
c_lon0, c_lon1 = (lon_m, lon1) if sj else (lon0, lon_m)
children.append(build_quadtree(q, c_lat0, c_lat1, c_lon0, c_lon1,
k_min, max_depth, depth + 1))
return {"bbox": (lat0, lat1, lon0, lon1), "children": children, "leaf": False}
# 3 通りの k_min で分割
for k_min in [2, 5, 10]:
root = build_quadtree(pts_all, *BBOX, k_min=k_min, max_depth=14)
leaves = [lf for lf in collect_leaves(root) if lf["pts"]]
print(f"k_min={k_min}: leaf={len(leaves)}")
''') + """
図と読み取り
なぜこの図か: セル境界の矩形と 点の分布を1枚に重ねることで、
「点が密な地域は細かい矩形、疎な地域は粗い矩形」という適応分割の本質が一目で分かる。
3つの k_min を並べると、k_min が大きいほど leaf が 合体して粗くなる過程が見える。
""" + figure("assets/L10_quadtree.png",
"図3: adaptive quadtree (k_min=2/5/10)。色は leaf 内点数。都市は細かく、山間は粗いセルに自動分割") + """
この図から読み取れること:
- 広島市・福山市付近は細かい矩形が密集: カメラが多いため細かく割れる
- 北部山間・島嶼部は大きな矩形: 点が疎なので粗いまま leaf 化
- k_min=2 → 5 → 10 と増やすと leaf 数が減る(矩形が大きくなる): 表4 の数値で裏付け
- すべての leaf が k_min 以上を保証: 矩形最小色 (薄黄) でも k_min 人いる
- 仮説H3 を支持: 固定グリッドが原理的に解けない「都市と山間の不均衡」を quadtree が解消
表と読み取り
表4: k_min ごとの quadtree 統計
""" + qt_summary_html + """
読み取り:
k_min を上げる → leaf 数が減る → 中央 leaf 面積が大きくなる という関係が定量化されている。
min 点数列は必ず k_min 以上で、設計通り 「全 leaf が k_min を満たす」保証が効いている。
分析5 のトレードオフ図で、これらの数値が 固定グリッドより左上(=同じ k を小さい平均面積で達成) に
配置されることを確認する。
"""),
# =========================================================================
# 7. 分析4: 差分プライバシ
# =========================================================================
("分析4: 差分プライバシ(Laplace ノイズ追加)", """
狙い
「全ての点に小さなランダムノイズを足したら、各点はどれくらい元位置からずれるか?」
を ε=1.0 → 0.5 → 0.1 と保護を強めながら可視化し、
「ε とジッタ距離は反比例する」という核心を体感する。仮説H4の検証。
ツールとしての差分プライバシ (DP):
- 入力: 351点の (lat, lon) と保護パラメータ ε(イプシロン)
- 出力: 351点の (lat+ノイズ, lon+ノイズ) — ノイズは Laplace 分布 から無作為抽出
- 保証 (黒箱で OK): 攻撃者が「あなた1人がデータに含まれていたかどうか」を見破る確率の比が e^ε 以下に抑えられる。
これが「数学的にプライバシを保証する」の意味で、内部の不等式 (Pr[出力∈S | あなた含む] / Pr[出力∈S | あなた含まない] ≤ e^ε) の証明は 気になる人向け。
- k-匿名性との根本的違い:
k-匿名性は 「同じセルに k 人いれば守られる」というグループ化保護で、攻撃者が背景知識を持つと壊れる。
DP は 「1人を入れ替えても出力分布がほぼ同じ」という確率的保護で、背景知識に強い。
- 限界: ε を強くすると(例 0.1) ジッタが大きすぎて地図が 使い物にならない。
また、同じデータに何回もクエリを投げると ε が累積消費される(composition theorem)。
手法 (Laplace 機構)
- 入力: 351点 + 感度 Δ(本レッスンでは 1km 相当) + ε ∈ {1.0, 0.5, 0.1}
- 処理:
- Laplace 分布のスケール
b = Δ / ε を決める (ε 小→b 大→ノイズ大)
- 緯度・経度それぞれに 独立に
Laplace(0, b) から抽出したノイズを加算
- ジッタ距離 (km) = √((Δlat·LAT_KM)² + (Δlon·LON_KM)²) で評価
- 出力: 各点のジッタ後位置 + 中央ジッタ・95%点・最大値の集計表
- パラメータの意味:
ε=1.0 → b=1km → 中央ジッタ ≈1.3km(緩い保護, 地図実用OK)。
ε=0.1 → b=10km → 中央ジッタ ≈13km(強保護, 県全域に拡散)
実装
""" + code('''
import numpy as np
# 感度 Δ = 「点を1km 動かす」を緯経度オフセットに換算
SENS_LAT = 1.0 / LAT_KM # 1km の緯度オフセット (≒0.0090度)
SENS_LON = 1.0 / LON_KM # 1km の経度オフセット (≒0.0109度)
rng = np.random.default_rng(42)
for eps in [1.0, 0.5, 0.1]:
b_lat = SENS_LAT / eps # Laplace スケール
b_lon = SENS_LON / eps
noise_lat = rng.laplace(0, b_lat, size=len(df)) # ノイズ抽出
noise_lon = rng.laplace(0, b_lon, size=len(df))
lat_dp = df["lat"].values + noise_lat # 加算 (核心1行)
lon_dp = df["lon"].values + noise_lon
# ジッタ距離 (km) を計算
d_km = np.hypot(noise_lat * LAT_KM, noise_lon * LON_KM)
print(f"eps={eps}: median={np.median(d_km):.2f}km, p95={np.percentile(d_km,95):.2f}km")
''') + """
ジャーゴン補足:
「Laplace 分布」=正規分布より 裾が長い(まれに大きなノイズが出る) 確率分布。
DP の Laplace 機構は数学的に「ε-DP を満たす最小ノイズ」として導かれる。
「再識別 (re-identification)」=匿名化したはずのデータから個人を特定し直すこと。
図と読み取り
なぜこの図か: オリジナルと 3つの ε を横一列に並べることで、
ε を強めるにつれて 赤点(DP後)がどれだけ 青点(オリジナル)から離れるかを直接比較できる。
赤線(対応線)はジッタの方向と距離を1点ずつ可視化する装置。
""" + figure("assets/L10_dp_noise.png",
"図4: Laplace ノイズによる差分プライバシ。ε=1.0(緩い)→0.5→0.1(強い)とジッタが10倍に拡大") + """
この図から読み取れること:
- ε=1.0: 中央ジッタ 1.3km、95%点 3.6km。「同じ市内にはずれるが県外には飛ばない」程度
- ε=0.5: 中央 2.8km、95%点 8.6km。市町境を越える点が増える
- ε=0.1: 中央 13km、95%点 38km。点が 県全域に拡散し、9% 弱は県外に飛ぶ。もはや「カメラ位置」と呼べる地理情報ではない
- 仮説H4 を支持: ε を 1.0→0.1 と 10倍強める と、ジッタも約10倍に拡大 (Laplace b = Δ/ε の単純な反比例)
表と読み取り
表5: ε ごとのジッタ統計
""" + dp_html + """
読み取り:
ε と Laplace b と中央ジッタは 1:1 対応の反比例。
ε=0.1 では「県内に収まる点」の割合さえ落ちる(最大ジッタが県外に飛ぶ)。
「ε はどう決めるか」は本来、データ提供者と利用者が 合意 すべき政策的な数字で、
教科書的には ε=1 前後が妥協値として使われることが多い。
"""),
# =========================================================================
# 8. 分析5: utility ↔ privacy トレードオフ
# =========================================================================
("分析5: utility ↔ privacy トレードオフ(全手法を1枚に)", """
狙い
「6種類の固定グリッド × 3種類の quadtree = 9手法」を同じ平面に置いて、
どの手法が utility-privacy のバランスで優位か を一目で比較する。
仮説H3 の最終確認 (quadtree が固定グリッドより左上に来るか)。
手法 (横軸=情報損失, 縦軸=プライバシ強度)
- 入力: 分析1-3 で計算済みの「セル面積中央 (km²)」と「セル内点数中央 (k)」
- 処理: x=セル面積(対数), y=k中央(対数) の散布。固定グリッドは丸、quadtree は星
- 結果の読み方:
左上=「面積小(=細かい)+k大(=保護強)」=理想。
右下=「面積大+k小」=最悪。
曲線がパレートフロンティア(達成可能境界)。
実装
""" + code('''
# 固定グリッド: 6 解像度
areas_km2 = [(r/1000)**2 for r in [100, 250, 500, 1000, 2000, 5000]]
median_k_per_res = [...] # 分析2 で計算済み
# quadtree: 3 通りの k_min
for k_min, _, leaves in QT_RESULTS:
leaf_areas = [(lf.lat1-lf.lat0)*LAT_KM * (lf.lon1-lf.lon0)*LON_KM
for lf in leaves]
leaf_sizes = [len(lf.pts) for lf in leaves]
med_area = np.median(leaf_areas) # 中央 leaf 面積
med_k = np.median(leaf_sizes) # 中央 k
# 同じ平面に重ね描画
ax.scatter(areas_km2, median_k_per_res, marker="o") # 固定グリッド
ax.scatter([med_area], [med_k], marker="*", c="red") # quadtree
ax.set_xscale("log"); ax.set_yscale("log")
''') + """
図と読み取り
なぜこの図か: utility(横軸) と privacy(縦軸) を 1枚の散布図に重ねることで、
「同じ k を達成するのに、どの手法が小さい面積で済むか」=パレート効率の比較ができる。
両軸とも対数にして広いダイナミックレンジを見やすく。
""" + figure("assets/L10_tradeoff.png",
"図5: 情報損失(セル面積) × プライバシ保護(k 中央) のトレードオフ。★=quadtree が同じ k に対して小さい面積") + """
この図から読み取れること:
- 固定グリッド(青〜黄丸)は y≈1 の水平線上: どの解像度でも k 中央が 1(351点の規模では占有セルがほぼ1点ずつ)。x軸方向(面積)だけが違う
- quadtree(★赤)は右側で y が上昇: 面積は大きいが、k 中央が 12〜76 と固定グリッドより圧倒的に大きい
- k_min=2 → 5 → 10 と上げると、★も右上に動く: k_min を厳しくすると leaf が合体して大きくなる
- 仮説H3 は部分支持: 同じ k を達成する固定グリッド設定がそもそも存在しない (固定グリッドはどこまで粗くしても k 中央=1) ので、quadtree の優位は「同じ k を達成できること自体」。
ただし、quadtree の中央面積は km単位の大きな値で、 「絶対的に細かい」わけではない。 351 点というサンプルサイズの少なさが構造的な制約を生んでいる
- 差分プライバシは別軸: DP は「セル面積」概念がない (=連続空間でジッタ) ので、この図には載らない。utility-privacy の構造が k-匿名性系とは異なる
表と読み取り
表6: 9手法の中央セル面積×k中央 一覧
""" + tradeoff_html + """
読み取り:
たとえば「k≧5 を確保したい」場合、固定グリッドでは 5km 格子でも k 中央=1 (=多くの点が孤立) で達成不能。
一方 quadtree (k_min=5) は k 中央=12 で確実に達成。
ただし quadtree の中央 leaf 面積は 898 km²(山間離島では1つの leaf が県全体の1/16級まで成長) で、
「点が疎な地域では適応分割でも結局粗くせざるをえない」という根本的な制約も同時に見える。
このトレードオフ表は 公開ガイドライン設計の数値根拠 になる
(例: 「県全域オープンデータでは k_min=5 を要求 → 山間部は 30km 級のセル面積を許容する」など)。
"""),
# =========================================================================
# 9. 仮説検証と考察
# =========================================================================
("仮説検証と考察", """
仮説と結果の照合
| # | 仮説 | 判定 | 根拠 |
|---|
| H1 |
固定グリッドは万能ではない (均一セル≠均一保護) |
支持 |
図1 で 1km格子でも山間部はほぼ k=1。表1 の k=1セル割合列が 250m で 70%超を示す。
均一格子では都市過剰粗・山間過剰細という構造的歪みが避けられない。 |
| H2 |
k を上げるほど大きなセルが必要 (k≧10 は km級) |
支持(かつ予想以上に厳しい) |
図2/表3 で k≧2 達成率は 1km で22%, 5km でも 67%にとどまる。
k≧5 は 5km 格子でも 24%、k≧10 は 2.8%。
曲線は単調に右上シフトし、351点規模では固定グリッドだけで k≧5 を多くの点で確保するのは原理的に困難。
これが quadtree (適応分割) が必要な強い動機。 |
| H3 |
quadtree は固定グリッドより utility-privacy で優位 |
部分支持 |
図5/表6 で★(quadtree)は k 中央 12〜76を達成、固定グリッドはどの解像度でも k 中央=1。
「k を確保できるかどうか」では quadtree の圧勝。
ただし quadtree の中央 leaf 面積は 200〜3500 km²と非常に大きく、
理論的に期待した「都市部の小セル」は表面化していない (中央値計算上、山間離島の巨大 leaf が支配)。
351 点というサンプル規模では適応分割でも限界があることを実証 — これは教育的な「失敗から学ぶ」点。 |
| H4 |
差分プライバシはε で揺らぎが反比例的に変わる |
支持 |
表5 で ε=1.0→0.5→0.1 と10倍強めると、中央ジッタも約10倍 (0.7km→1.4km→7km)。
Laplace b=Δ/ε の数学的関係が実データでも素直に現れた(中央ジッタ 1.3km→2.8km→13km)。
ε=0.1 では ジッタが県全域級で「地図用途には使えない」も確認。 |
| H5 |
1棟特定は10m格子で起きる (GPS 精度と一致) |
支持 |
下記 補助表7 で 10m 格子では占有セルの 95%超が k=1。
GPSの生精度と符合し、「位置データは本質的に強い識別子」を数値で確認できた。 |
補助: 家1棟スケールでの k=1 リスク (仮説H5の根拠)
""" + house_html + """
5〜10m 格子では占有セルの 90% 超が k=1 (=同一セル内に他点なし)。
これは 1棟単位での特定が技術的に可能 な粒度。
災害時に「自宅前の道路カメラ」と一意に紐づく公開は、住民属性データと重なれば
プライバシ侵害につながりうる。
考察
- 「均一格子の盲点」と quadtree の意義:
H1 の支持が示すように、固定グリッドは「均一」と「公平」を取り違える典型。
「均一に粗く・均一に細かく」では均一に保護できない。
quadtree のような適応的データ依存分割が、現代の位置情報公開では実質標準になりつつある。
- k-匿名性 vs 差分プライバシ — 哲学が異なる:
k-匿名性は「グループ化による保護」で、攻撃者が背景知識を持つと壊れる(homogeneity / background-knowledge attack)。
差分プライバシは「数学的な保証」: 1点を入れ替えても出力分布の比が e^ε に収まる。
ただし ε=1.0 でも数km単位のジッタが出るので、地図公開とは相性が悪い場面もある。
実務では ハイブリッド(quadtree + DP) が増えている。
- 「災害時の例外」:
図4 ε=0.1 のように大きく揺らした地図は、避難所への誘導には使い物にならない。
緊急時の精度と 平時のプライバシは両立しないこともある。
これは 政策論(緊急時例外規定など)で扱われる領域。
データ屋は選択肢のコストを提示する役で、最終判断は政策側。
- "1棟特定" 閾値の心得:
表7 から 10m 格子で 95% 以上が k=1。
GPS の生精度が 5〜10m であることと符合する。
「位置データは本質的に強い識別子」を、自分の手で数字として確認することが、
心得カテゴリの中核。
- 方法の限界:
本レッスンの感度 Δ=1km は「1km 動かす」という人為的な決め方。
実務では「世帯1軒の入れ替え」「個人1人の入れ替え」を 感度と定義する。
そこの設計が ε と並んで DP 適用の成否を決める鍵。
"""),
# =========================================================================
# 10. 発展課題
# =========================================================================
("発展課題(結果から導かれる新たな問い)", """
各課題は、上の 結果 と 新たな仮説 に裏打ちされている。
「結果X→新仮説Y→課題Z」の3段で記述。
- l-diversity / t-closeness — k-匿名性の弱点を補う
- 結果X: 分析1/3 でk-匿名性は「人数」だけを見るが、セル内属性が 均質(全員が同じ管理区分=道路カメラ etc.) なら k≧5 でも事実上特定できる(homogeneity attack)
- 新仮説Y:
camera_list.csv の 路河川名や 管理区分(道路/河川) を属性として、各セルに l 種類以上の属性を要求すれば、homogeneity に強くなる
- 課題Z: 各セルに「ユニーク管理区分数 ≥ 2」の制約を追加した quadtree を実装。 leaf 数と k 中央がどう変わるかを比較
- quadtree の k_min 自動調整 — データ依存決定
- 結果X: 分析3 で k_min を 2/5/10 と固定値で振ったが、最適値はデータの密度分布に依存する
- 新仮説Y: 全点数 N に対して k_min ≈ √N(本データなら ≈19) が utility-privacy バランスの理論的最適に近いはず
- 課題Z: k_min=√N で再実行し、図5 の★位置が他の k_min より 左上に来るかを検証
- 差分プライバシの累積予算 — composition theorem の可視化
- 結果X: 分析4 は1回のクエリだけを扱った。実務では同じデータに 複数回クエリを投げる
- 新仮説Y: 同じデータセットに5回クエリを投げると ε が 5倍に消費(strong composition なら √5 倍)。元 ε=1.0 でも 累積 ε=5 となり実質保護が弱い
- 課題Z: 同じノイズ生成を5回繰り返してその度の ε 消費を可視化、「クエリ予算」の概念を体験する
- quadtree + DP のハイブリッド — 2段保護の実装
- 結果X: 分析3/4 はそれぞれ単独で長所短所がある。 quadtree は背景知識攻撃に弱く、DP は精度を犠牲にしすぎる
- 新仮説Y: 「まず quadtree でセル化 → 各セルの 件数 に Laplace ノイズを加える」2段保護なら、両方の長所を取れるはず
- 課題Z: quadtree の各 leaf 件数に
Laplace(0, 1/ε) を加算した「ぼかし件数地図」を作る。 utility(=件数誤差) と privacy(ε) のトレードオフを再描画
- 政策接続 — 自治体オープンデータの公開ガイドライン調査
- 結果X: 表7 から「10m 格子では 95% が k=1」=技術的に1棟特定可能。 これがどの法規でカバーされるかは未確認
- 新仮説Y: 個人情報保護法・統計法・各自治体条例のいずれかに、「位置データの最低粒度」を定める条文があるはず
- 課題Z: 広島県・国の公開ガイドラインを読み、「位置情報を含むオープンデータの匿名化基準」がどう定められているかを調査・引用付きでまとめる
"""),
]
html = render_lesson(
num=10,
title="プライバシ・グリッド — k-匿名性 / quadtree / 差分プライバシ",
tags=["v2-rewrite", "心得", "応用基礎", "プライバシ", "ガバナンス"],
time="120分",
level="心得 / 応用基礎",
data_label='DoBoX #1279 '
'カメラ位置 351点 (位置情報のプライバシ保護シミュレーション)',
sections=sections,
script_filename="lessons/L10_privacy_grid.py",
)
out = LESSONS / "L10_privacy_grid.html"
out.write_text(html, encoding="utf-8")
print(f"\nsaved: {out}")
print(f" PNGs: 5 generated in {ASSETS}")
for p in sorted(ASSETS.glob("L10_*.png")):
print(f" - {p.name} ({p.stat().st_size//1024} KB)")
print(f" CSVs: generated in {ASSETS}")
for p in sorted(ASSETS.glob("L10_*.csv")):
print(f" - {p.name} ({p.stat().st_size//1024} KB)")
