導入
センサや観測点がまばらなとき、よくやりたくなるのが「地図にする」ことです。
- 環境（気温・PM2.5・騒音・水質）
- 医工学（病院内の混雑・機器温度・電力）
- 実験（反応条件の空間分布、装置内の温度ムラ）
- 都市・施設（混雑、遅延、通信品質）
でも、観測点が少ないのに **“滑らかな地図”**だけを出すと危険です。
なぜなら、地図のどこが「ちゃんと観測に支えられた推定」で、どこが「だいぶ当てずっぽう」かが見えないからです。
そこで本記事では、
- 空間補間（Spatial interpolation）で推定地図を作り、
- さらに 不確実性（Uncertainty） を地図として出して、
- 「推定値」と「自信のなさ」をセットで提示する
最小実験を Google Colab で体験します。
このデモは教育用で、地球物理や大気モデルなどの厳密再現ではありません。
ただし「空間補間＋不確実性の可視化」という方法論は、多くの現場でそのまま役立ちます。
TL;DR
- 空間補間は「見た目の地図」を作れるが、どこが信用できるかが重要。
- IDW（Inverse Distance Weighting）は簡単で速いが、不確実性は出しにくい。
- Kriging ≒ ガウス過程回帰（Gaussian Process Regression）は、
予測の 平均（推定値） と 標準偏差（不確実性） を同時に出せる。 - 不確実性を「別の地図」にするだけでなく、推定地図の透明度に反映すると“映える”し、誤解も減る。
1. 今日やること（最小構成）
このデモでは、2次元平面（1km四方のイメージ）に「真の場（ground truth）」を作り、ランダムな観測点だけを与えます。
- 真の場（見えないはずの分布）を合成
- まばらな観測点（センサ）を置き、ノイズ付き観測を得る
- IDWで補間（ベースライン）
- ガウス過程（Kriging）で補間し、不確実性も地図化
- “地図が映える”見せ方（不確実性を透明度へ）
2. 用語を学部生向けに（重要）
2.1 空間補間（Spatial interpolation）
観測点がまばらでも、全域の値を推定してグリッド（格子）上の地図にします。
「どこでも値がある地図」は直感的ですが、観測が無い場所ほど危ういです。
2.2 不確実性（Uncertainty）
ここで言う不確実性は「誤差そのもの」ではなく、
観測点の配置とモデル仮定（滑らかさ等）から見て、
どれくらい自信を持って推定できるか
の指標です。
2.3 Kriging と ガウス過程（GP）の関係
Kriging は地質統計（geostatistics）で有名な空間補間手法です。
実は「ある種の仮定のもとで、ガウス過程回帰とほぼ同じ形」になります。
本記事では実装しやすい scikit-learn の GaussianProcessRegressor を使い、
- 予測平均（mean）＝推定値の地図
- 予測標準偏差（std）＝不確実性の地図
を作ります。
3. Google Colabで動かす（コピペでOK）
以下を Colab に上から貼り付けて実行してください。
fig/ に図を保存します。
セル1：準備
import os
import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel
from sklearn.preprocessing import StandardScaler
SEED = 42
rng = np.random.default_rng(SEED)
os.makedirs("fig", exist_ok=True)
セル2：真の場（ground truth）を合成し、観測点を作る
def build_true_field():
"""滑らかな2D場（山が2つ＋うねり）を合成。※物理モデルではないデモ用"""
def f(x, y):
# 2つのガウス山 + ゆるい波
g1 = 1.3*np.exp(-(((x-0.75)/0.12)**2 + ((y-0.35)/0.10)**2))
g2 = 0.9*np.exp(-(((x-0.25)/0.18)**2 + ((y-0.75)/0.12)**2))
w = 0.15*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)
return g1 + g2 + w
return f
f_true = build_true_field()
# グリッド（地図を描く範囲）
N_GRID = 160
xs = np.linspace(0, 1, N_GRID)
ys = np.linspace(0, 1, N_GRID)
XX, YY = np.meshgrid(xs, ys)
Z_true = f_true(XX, YY)
# 見やすいように 0〜100 にスケーリング
mn, mx = Z_true.min(), Z_true.max()
Z_true_01 = (Z_true - mn) / (mx - mn)
Z_true_100 = 100 * Z_true_01
# 観測点（センサ）をランダム配置
N_SENSORS = 60
X_obs = rng.uniform(0, 1, size=(N_SENSORS, 2)) # (x,y)
z_clean = 100 * (f_true(X_obs[:,0], X_obs[:,1]) - mn) / (mx - mn)
# 観測ノイズ（測定誤差のイメージ）
NOISE_SIGMA = 4.0
z_obs = z_clean + rng.normal(0, NOISE_SIGMA, size=N_SENSORS)
print("z_obs range:", float(z_obs.min()), float(z_obs.max()))
セル3：図1（真の場＋観測点）※デモ用の“答え合わせ”
plt.figure(figsize=(7.5,6))
cf = plt.contourf(XX, YY, Z_true_100, levels=30)
plt.colorbar(cf, label="value (0-100, synthetic)")
plt.scatter(X_obs[:,0], X_obs[:,1], s=35, edgecolor="k", linewidth=0.5, label="sensors")
plt.title("Ground truth (toy demo) + sensor locations")
plt.xlabel("x"); plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig1_ground_truth_and_sensors.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved: fig/fig1_ground_truth_and_sensors.png")
セル4：ベースライン（IDW補間：Inverse Distance Weighting）
「近い観測点ほど重く、遠いほど軽く」する単純な補間です。
不確実性は出しにくいですが、最初の比較として便利です。
def idw_predict(X_obs, z_obs, X_pred, power=2.0, k=12, eps=1e-12):
"""IDW: 予測点X_predごとに、近いk点の距離の逆数^powerで重み付け平均"""
# 距離行列: (n_pred, n_obs)
d = np.sqrt(((X_pred[:,None,:] - X_obs[None,:,:])**2).sum(axis=2)) + eps
# k近傍だけ使う
idx = np.argpartition(d, kth=k-1, axis=1)[:, :k]
d_k = np.take_along_axis(d, idx, axis=1)
z_k = z_obs[idx]
w = 1.0 / (d_k**power)
w = w / (w.sum(axis=1, keepdims=True) + eps)
return (w * z_k).sum(axis=1)
# 予測点（グリッド）
X_grid = np.column_stack([XX.ravel(), YY.ravel()])
Z_idw = idw_predict(X_obs, z_obs, X_grid, power=2.0, k=12).reshape(XX.shape)
plt.figure(figsize=(7.5,6))
cf = plt.contourf(XX, YY, Z_idw, levels=30)
plt.colorbar(cf, label="IDW predicted value")
plt.scatter(X_obs[:,0], X_obs[:,1], s=25, edgecolor="k", linewidth=0.4)
plt.title("IDW interpolation (baseline)")
plt.xlabel("x"); plt.ylabel("y")
plt.tight_layout()
plt.savefig("fig/fig2_idw.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved: fig/fig2_idw.png")
セル5：ガウス過程（Kriging）で平均と不確実性を出す
ポイント：
- Maternカーネル：滑らかすぎない現実寄りの相関を表現しやすい
- WhiteKernel：測定ノイズ（“ナゲット”に相当）の表現
-
return_std=Trueで予測標準偏差を得る
# 入力スケールの安定化（座標を標準化）
scaler = StandardScaler()
X_obs_s = scaler.fit_transform(X_obs)
X_grid_s = scaler.transform(X_grid)
kernel = ConstantKernel(1.0, (1e-2, 1e2)) * Matern(
length_scale=[1.0, 1.0],
length_scale_bounds=(1e-2, 1e2),
nu=1.5
) + WhiteKernel(noise_level=1.0, noise_level_bounds=(1e-5, 1e2))
gp = GaussianProcessRegressor(
kernel=kernel,
normalize_y=True,
n_restarts_optimizer=3,
random_state=SEED
)
gp.fit(X_obs_s, z_obs)
mu, std = gp.predict(X_grid_s, return_std=True)
Z_gp_mean = mu.reshape(XX.shape)
Z_gp_std = std.reshape(XX.shape)
print("Learned kernel:", gp.kernel_)
print("mean range:", float(Z_gp_mean.min()), float(Z_gp_mean.max()))
print("std range:", float(Z_gp_std.min()), float(Z_gp_std.max()))
セル6：図3（GPの推定値の地図）
plt.figure(figsize=(7.5,6))
cf = plt.contourf(XX, YY, Z_gp_mean, levels=30)
plt.colorbar(cf, label="GP mean (predicted value)")
plt.scatter(X_obs[:,0], X_obs[:,1], s=25, edgecolor="k", linewidth=0.4)
plt.title("Gaussian Process / Kriging: mean (estimate)")
plt.xlabel("x"); plt.ylabel("y")
plt.tight_layout()
plt.savefig("fig/fig3_gp_mean.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved: fig/fig3_gp_mean.png")
セル7：図4（GPの不確実性の地図：標準偏差）
plt.figure(figsize=(7.5,6))
cf = plt.contourf(XX, YY, Z_gp_std, levels=30)
plt.colorbar(cf, label="GP std (uncertainty)")
plt.scatter(X_obs[:,0], X_obs[:,1], s=25, edgecolor="k", linewidth=0.4)
plt.title("Gaussian Process / Kriging: std (uncertainty map)")
plt.xlabel("x"); plt.ylabel("y")
plt.tight_layout()
plt.savefig("fig/fig4_gp_std.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved: fig/fig4_gp_std.png")
セル8：図5（“映える”見せ方：不確実性を透明度にする）
「推定値の色」は同じでも、不確実性が大きい領域を薄くすると、
“分かっている場所／分かっていない場所”が直感的になります。
# 不確実性を 0〜1 に正規化して、alpha（透明度）に変換
std_norm = (Z_gp_std - Z_gp_std.min()) / (Z_gp_std.max() - Z_gp_std.min() + 1e-12)
alpha = 1.0 - 0.75*std_norm # 不確実性が大きいほど薄く
plt.figure(figsize=(7.5,6))
im = plt.imshow(
Z_gp_mean,
origin="lower",
extent=[0,1,0,1],
alpha=alpha
)
plt.colorbar(im, label="GP mean (predicted value)")
plt.scatter(X_obs[:,0], X_obs[:,1], s=25, edgecolor="k", linewidth=0.4)
plt.title("Mean map with uncertainty as transparency (visually safer)")
plt.xlabel("x"); plt.ylabel("y")
plt.tight_layout()
plt.savefig("fig/fig5_mean_with_uncertainty_alpha.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved: fig/fig5_mean_with_uncertainty_alpha.png")
セル9：図6「次にセンサを置くならここ」（不確実性上位点）
これは “空間補間＋不確実性” の強い使いどころです。
不確実性が高い場所を優先して観測すると、地図の品質が上がりやすいです。
# 不確実性が高い上位K点を提案
K = 12
flat_idx = np.argsort(Z_gp_std.ravel())[::-1][:K]
X_next = X_grid[flat_idx]
plt.figure(figsize=(7.5,6))
cf = plt.contourf(XX, YY, Z_gp_std, levels=30)
plt.colorbar(cf, label="GP std (uncertainty)")
plt.scatter(X_obs[:,0], X_obs[:,1], s=25, edgecolor="k", linewidth=0.4, label="current sensors")
plt.scatter(X_next[:,0], X_next[:,1], s=80, marker="*", edgecolor="k", linewidth=0.8, label="suggested next")
plt.title("Where to sample next? (top uncertainty points)")
plt.xlabel("x"); plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig6_next_sampling.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved: fig/fig6_next_sampling.png")
セル10：図をまとめてダウンロード（Colab用）
import shutil
shutil.make_archive("figures", "zip", "fig")
print("Created figures.zip")
4. 実行結果の図
図1：真の場＋観測点（デモ用の答え合わせ）
図2：IDW（ベースライン補間）
図3：GP/Kriging 推定値（平均）
図4：GP/Kriging 不確実性（標準偏差）
図5：“映える”安全な地図（不確実性を透明度へ）
図6：次に観測すべき候補（不確実性上位）
5. 図の読み方（推定値と不確実性をセットで読む）
5.1 推定値（図3）は「地図」だが、単体では危険
図3だけ見ると、どこでもそれっぽい値が塗られていて安心してしまいます。
しかし実際は、観測点が少ない場所では「そう見えるだけ」の可能性があります。
5.2 不確実性（図4）は「どこが分かっていないか」の地図
一般に、不確実性が高くなりやすいのは：
- 観測点から遠い場所
- 境界付近（外挿に近い）
- 観測点配置が偏っている領域
です。
図4が明るい場所は「推定値が間違っている」ではなく、
「推定値に自信がない」という意味です。
5.3 図5（透明度）は、誤解を減らす“伝え方”
地図を使うと「色＝事実」に見えやすいです。
そこで、不確実性が大きい場所を薄くすることで、
- “ここは確度が高い”
- “ここは推測が強い”
が一目で伝わり、意思決定が安全になります。
6. 実務で使うときの注意点（ここがインフォマティクス）
6.1 不確実性は“モデル仮定の上”にある
GP/Krigingの不確実性は、主に次に依存します。
- 相関の仮定（カーネル、長さスケール）
- 定常性・等方性（場所によらず同じ性質とみなす等）
- ノイズ（測定誤差）
仮定が外れると、不確実性が過小評価・過大評価されることがあります。
「不確実性が小さいから絶対に正しい」ではなく、“仮定の中で自信がある”です。
6.2 座標が緯度経度なら距離の扱いに注意
本記事は 0〜1 の平面座標として扱いました。
緯度経度（lat/lon）をそのままユークリッド距離で扱うと歪む場合があります。
狭い範囲なら近似できても、広域なら投影や測地距離を検討します。
6.3 外部特徴（標高、土地利用、風向など）があると強い
純粋な空間相関だけでなく、説明変数を加えたモデル（回帰＋GP残差など）にすると、
「なぜその分布になるか」も含めて推定が安定します。
まとめ
本記事では、ダミーデータを使って
- 空間補間（IDW / GP=Kriging）
- 不確実性（GP std）
- “映える”かつ誤解を減らす可視化（透明度）
を Google Colab で再現しました。
「地図を作る」だけで終わらず、
地図の信用度（不確実性）を一緒に出すと、研究でも運用でも意思決定が強くなります。
次は、あなたのデータに合わせて
- 観測点の偏り
- ノイズの大きさ
- 長さスケール（どれくらい滑らかか）
- 次にどこを測るべきか
を調整し、現場にフィットした“見せ方”を作ってみてください。