位置情報の誤差を含むガウス過程モデルにおける事後推論：ウォーカー湖のウラニウムデータを用いた事例分析

背景と研究目的

鉱業分野において空間確率モデルは強力なツールとなります。

• 用途: 鉱資源探査において、ボーリングによる地質試料の分析から貴金属鉱石の有無や濃度を推測します。
• 課題: データに強い空間相関があるため、観測できる地下状態のごく一部のみです。これにより完全な地球物理モデルの構築は困難になります。
• 解決策: 近傍のデータを用いて、新しい地点における関心変数を予測する確率モデルを構築します。

モデルとデータの概要

この事例では、ガウス過程モデルが位置誤差を含む場合の扱い方を示します。

データソース

• ウォーカー湖（Walker Lake）のウラニウム・バナジウムデータ
  • 点ごとの濃度測定値を使用
  • 出典: Isaaks & Srivastava の『応用地統計学の入門』
  • パッケージ:

    gstat

    に付属

モデル設定

位置情報の誤差は共分散行列を変化させ、推定問題そのものを変える要因となります。

• 記号定義:
  • (\tilde{\mathbf{s}}_i): 記録された座標（観測地点）
  • (\mathbf{s}_i): 実際には観測が起きた潜在座標
  • (\Delta_i): 位置誤差 ((\Delta_i \sim \operatorname{Normal}(\mathbf{0}, \sigma_s^2 I_2)))
• 評価点: ガウス過程を記録された座標 (\tilde{\mathbf{s}}_i) の代わりに、実際の座標 (\mathbf{s}_i = \tilde{\mathbf{s}}_i + \Delta_i) で評価します。
• 事前分布の仮定: スケールパラメータ (\sigma_s) を既知量として扱い、異なる誤差レベルをモデルに反映させています。

確率モデル（ベイズモデル）

以下が導入されたモデルプロセスです：

\begin{aligned} 
\Delta_i &\sim \mathrm{Normal}(\mathbf{0}, \sigma_s^2\mathbf{I}_2) \\ 
\mathbf{s}_i &= \tilde{\mathbf{s}}_i + \Delta_i \\ 
\mu &\sim \mathrm{Normal}(2, 2) \\ 
\sigma &\sim \mathrm{HalfNormal}(1) \\ 
\ell &\sim \mathrm{HalfNormal}(100) \\ 
\sigma_0 &\sim \mathrm{HalfNormal}(0.5) \\ 
f(\cdot) \mid \sigma,\ell &\sim \mathcal{GP}(0, \sigma^2 c(\cdot,\cdot;\ell)) \\ 
Y_i \mid f,\mu,\sigma_0,\mathbf{s}_i &\sim \mathrm{Normal}(\mu + f(\mathbf{s}_i), \sigma_0) 
\end{aligned}

特徴:

• 固定位置のモデルと比較して、潜在座標が変化するため共分散行列も変化します（計算上困難）。
• 観測点における潜在ガウス過程の値を積分化（マージナライズ）するために、

  pm.gp.Marginal

  を使用します。

実装とシミュレーション実験

ノイズレベルを増加させたデータセットを作成し、モデルのパラメータ推定量がどのように変化するかを検査します。

1. データの前処理と摂動生成

元の座標に増加するノイズを加えて摂動を与えます（

multipliers = [12.0, 25.0, 40.0]

from pathlib import Path
import arviz as az
import matplotlib as mpl
import numpy as np
import pandas as pd
import pymc as pm
# ... (他のインポート)

RANDOM_SEED = 8927
np.random.seed(RANDOM_SEED)

# データ読み込み
df = pd.read_csv(Path("../../data/walker/cleaned.csv"))
X_walker = df[["x_m", "y_m"]].to_numpy()
y_walker = df["u_log10_p1"].to_numpy()

# 予測グリッドの作成
n_prediction_grid = 70
buffer_fraction = 0.1
x_range = df["x_m"].max() - df["x_m"].min()
x_limits = (df["x_m"].min() - buffer_fraction * x_range, df["x_m"].max() + buffer_fraction * x_range)
# ... (グリッド生成コード)

# 摂動データの生成
eps = np.random.randn(*X_walker.shape)
multipliers = [12.0, 25.0, 40.0]
noisy_xs = [X_walker + m * eps for m in multipliers]

2. モデルの構築 (

PyMC

)

pm.Data

location_coords = {"obs": np.arange(X_walker.shape[0]), "coord": ["x", "y"]}
selected_location_idx = np.sort(np.random.default_rng(RANDOM_SEED).choice(X_walker.shape[0], size=4, replace=False))
location_error_idatas = {}

with pm.Model(coords=location_coords) as location_error_gp_model:
    X_noisy = pm.Data("X_noisy", noisy_xs[0], dims=("obs", "coord"))
    y = pm.Data("y", y_walker, dims="obs")
    σ_s = pm.Data("σ_s", multipliers[0])

    # 潜在座標の定義
    Δs = pm.Normal("Δs", 0.0, σ_s, dims=("obs", "coord"))
    X_true = pm.Deterministic("X_true", X_noisy + Δs, dims=("obs", "coord"))

    # 事前分布
    μ = pm.Normal("μ", 2.0, 2.0)
    σ = pm.HalfNormal("σ", 1.0)
    ℓ = pm.HalfNormal("ℓ", 100.0)
    σ0 = pm.HalfNormal("σ0", 0.5)

    # ガウス過程定義
    cov = σ**2 * pm.gp.cov.Matern52(2, ls=ℓ)
    gp_location = pm.gp.Marginal(mean_func=pm.gp.mean.Constant(μ), cov_func=cov)
    
    # 事後確率の計算
    gp_location.marginal_likelihood("y_obs", X=X_true, y=y, sigma=σ0, dims="obs")

    # シミュレーションループ
    for multiplier, X_noisy_value in zip(multipliers, noisy_xs):
        pm.set_data({"X_noisy": X_noisy_value, "σ_s": multiplier})
        location_error_idatas[multiplier] = pm.sample(
            chains=2,
            cores=2,
            target_accept=0.95,
            mp_ctx="spawn",
            random_seed=RANDOM_SEED,
            progressbar=False,
        )

3. サンプリング結果の診断

高いノイズレベルに対して「分岐（divergences）」や「チェーン深度の限界」に関する警告が表示される場合があります。

注：分岐回数の増加は、高次元な後方分布や複雑な共分散構造による探索の困難さを示唆しています。

4. 事後予測と可視化

予測グリッドに対する事後予測を計算し、ベイズガウス過程モデルとナイーブな核密度推定（KDE）を比較します。

# ... (前処理略)

location_surface_grids = {}
with location_error_gp_model:
    for multiplier, X_noisy_value in zip(multipliers, noisy_xs):
        pm.set_data({"X_noisy": X_noisy_value, "σ_s": multiplier})
        
        # 事後平均の抽出
        posterior_mean_point = {
            name: location_error_idatas[multiplier].posterior[name].mean(("chain", "draw")).values
            for name in ["μ", "σ", "ℓ", "σ0", "Δs", "X_true"]
        }
        
        # グリッド上の予測
        f_mean, _ = gp_location.predict(Xnew, point=posterior_mean_point, diag=True, pred_noise=False)
        location_surface_grids[multiplier] = f_mean.reshape(n_prediction_grid, n_prediction_grid)

# ナイーブな KDE の計算 (比較用)
naive_kde_grids = {}
kde_bandwidths = {
    multiplier: location_error_idatas[multiplier].posterior["ℓ"].mean(("chain", "draw")).item()
    for multiplier in multipliers
}
# ... (KDE 計算と可視化コード略)

fig.savefig(figure_dir / "error-in-location-grid.png")

図 6：結果の解釈

ガウス過程モデルにおける点座標が不確実である場合のデータと事後推論

• 上段行:
  • 空心円：実在する元の座標
  • 着色された点：ノイズの混ざった観測データ（ウラニウム濃度に応じた色）
  • 大きな円：座標誤差の空間スケールを示す
• 2 段目:
  • 4 つの選ばれた点に対する事後高確率密度領域の輪郭線
  • (\sigma_s)（位置誤差）が増加するにつれて、これらの領域が滑らかに拡大していることが確認できます。

発見と結論

1. 特徴の保存: 座標の不確かさが広がるにもかかわらず、基礎となる表面の特徴の一部が保存されていることが観察されます。
   • 例：左下の明るい部分や左上・右下の暗い地域など。
   • 摂動の深刻さを考慮しても、こうしたデータで分析を進行させることは驚くべき成果です。
2. ナイーブなアプローチとの比較:
   • バンド幅を事後平均スケール長さに設定した Nadaraya-Watson smoother を使用した単純なアプローチと比較しました。
   • その結果、単純な平滑化手法は空間的な変動を表す能力が乏しく、単なる大まかな平均しか提供できないことが明らかになりました。

要約: ガウス過程モデルを適切に改変することで、位置情報の誤差を伴うデータでも信頼できる空間推論が可能であることを実証しました。