데이터 출처가 불명? 베이지안 모델링으로 미지 좌표 분석

출처: Hacker News

광산 산업에서 공간 확률 모델을 사용하는 강력한 동기 사례가 있다. 광물 자원을 탐사할 때, 탐사자는 구멍을 뚫어 시료를 채취하고 그 물질에 귀중한 광물이 존재하거나 농도가 높은지를 검사한다. 이러한 데이터는 일반적으로 강한 공간 상관성을 보이지만, 지하 상황을 거의 관측할 수 없기 때문에 완전한 지구물리학 모델을 구축하는 것이 때때로 불가능하다. 다행히도 지표 투과 레이더와 중력 측정과 같은 원격 탐사 기술의 등장으로 지구 내부를 특성화하는 능력이 크게 향상되었다. 이 문제를 해결하기 위해 우리는 인근 데이터를 이용해 새로운 위치에서 관심 변수의 값을 예측하는 확률 모델을 만들고자 한다.

문제를 더 명확히 설명하기 위해, Walker Lake에서 측정된 우라늄과 바나듐 농도에 대한 점 데이터셋을 사용한다. 이 데이터는 Isaaks와 Srivastava의 An Introduction to Applied Geostatistics에서 유래했으며 R 패키지 gstat와 함께 제공된다.

지금까지 로봇공학, 공간 통계, 신경과학 등에서 가우시안 프로세스 모델이 어떻게 사용되는지에 대한 많은 예시를 보았을 것이다. 이제는 실제 데이터 포인트의 위치가 정확히 알려져 있지 않고 상당한 측정 잡음으로만 관측되는 경우를 다루는 보다 이색적인 예시를 살펴본다. 위치 오류는 공분산과 예측 문제 자체를 변화시키며, 이는 위치 오류에 관한 지구통계학 연구와 잡음이 섞인 공간 입력에 대한 GP 회귀 연구에서 강조된 바 있다(Cressie와 Kornak, Cervone와 Pillai).

처음에는 다소 특이하게 보일 수 있지만, 이 예시는 적절한 사전분포를 갖춘 베이지안 모델링이 모델의 거의 모든 부분을 가정에 맞게 수정·변경할 수 있음을 보여준다. 그런 다음 몬테카를로 방법을 사용해 추론을 수행하고 신뢰할 수 있는 파라미터 추정값을 얻는다.

몇 가지 표기법을 추가하면, (\tilde{\mathbf{s}}_i)는 기록된 좌표, (\mathbf{s}_i)는 실제 측정이 일어난 잠재 좌표를 의미한다. 우리는 (\mathbf{s}_i = \tilde{\mathbf{s}}_i + \Delta_i) 로 두고, (\Delta_i \sim \operatorname{Normal}(\mathbf{0}, \sigma_s^2 I_2)) 로 가정한다. 그리고 가우시안 프로세스를 (\tilde{\mathbf{s}}_i)가 아니라 (\mathbf{s}_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} ]

잠재 좌표가 바뀔 때마다 공분산 행렬이 변하기 때문에, 이 모델은 고정 위치 GP보다 계산적으로 더 어렵다. 우리는 pm.gp.Marginal을 사용해 관측값에 대한 잠재 GP 값을 적분한다.

우리는 잡음 수준을 점점 높여가며 데이터셋을 만들고, 모델 파라미터 추정값이 어떻게 변하는지 살펴볼 것이다. 먼저 원래 좌표에 점진적으로 잡음을 추가한다:

from pathlib import Path

import arviz as az
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
from matplotlib import colors
from matplotlib.patches import Circle

from stsp import use_sepia_gunmetal

RANDOM_SEED = 8927
np.random.seed(RANDOM_SEED)

theme, sepia_cmap = use_sepia_gunmetal()

plot_bg_color = "#1c1c1d"
plot_text_color = "#e8e8e8"

plt.style.use("dark_background")
mpl.rcParams.update(
    {
        "figure.facecolor": plot_bg_color,
        "axes.facecolor": plot_bg_color,
        "savefig.facecolor": plot_bg_color,
        "savefig.edgecolor": plot_bg_color,
        "text.color": plot_text_color,
        "axes.labelcolor": plot_text_color,
        "xtick.color": plot_text_color,
        "ytick.color": plot_text_color,
        "axes.edgecolor": "#828282",
    }
)

theme["paper"] = plot_bg_color

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()
y_range = df["y_m"].max() - df["y_m"].min()
x_limits = (df["x_m"].min() - buffer_fraction * x_range, df["x_m"].max() + buffer_fraction * x_range)
y_limits = (df["y_m"].min() - buffer_fraction * y_range, df["y_m"].max() + buffer_fraction * y_range)
x_new = np.linspace(*x_limits, n_prediction_grid)
y_new = np.linspace(*y_limits, n_prediction_grid)
x_new_mesh, y_new_mesh = np.meshgrid(x_new, y_new)
Xnew = np.column_stack([x_new_mesh.ravel(), y_new_mesh.ravel()])

eps = np.random.randn(*X_walker.shape)
multipliers = [12.0, 25.0, 40.0]
noisy_xs = [X_walker + m * eps for m in multipliers]

잡음이 추가된 좌표는 Figure X에 표시된다. 이 분석에서는 목표 변수는 그대로 유지한다. 다음으로 pm.Data 객체를 컨테이너로 사용해 좌표를 쉽게 교체할 수 있도록 모델을 구성한다.

location_coords = {"obs": np.arange(X_walker.shape[0]), "coord": ["x", "y"]}
selec

ted_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,
        )

location_error_parameter_summaries = {
    multiplier: az.summary(idata, var_names=["μ", "σ", "ℓ", "σ0"])
    for multiplier, idata in location_error_idatas.items()
}
location_error_diagnostics = pd.DataFrame(