不知道您的数据来自哪里？未知坐标的贝叶斯建模

使用空间概率模型的一个特别强烈的动机来自采矿业。在矿产资源勘探过程中，探矿者将通过钻孔采集地质样本并检查所得材料是否存在或浓缩有价值的矿石。这些数据通常显示出很强的空间相关性，但构建全面详细的地球物理模型有时是不可行的，因为我们只能观察到很少的地下条件，尽管探地雷达和重力测量等遥感技术的出现极大地提高了我们描述地球地下特征的能力。为了应对这一挑战，我们希望构建一个概率模型，该模型使用附近的数据来预测新位置处的感兴趣变量。为了更好地说明问题，我们将使用沃克湖的铀和钒点参考浓度测量数据集。数据源自 Isaaks 和 Srivastava 的《应用地统计学简介》，并随 R 包 gstat 一起分发。到目前为止，您可能已经看到了很多关于如何在机器人、空间统计、神经科学等领域使用高斯过程模型的示例。现在，我们将讨论一个更奇特的示例，该示例修改高斯过程模型以适应数据点的实际位置未知且仅在大量测量噪声下观察到的情况。空间位置误差改变了协方差和预测问题本身，这是 Cressie 和 Kornak 以及 Cervone 和 Pillai 在位置误差的地统计工作和噪声空间输入的 GP 回归工作中强调的一点。乍一看，这似乎是一个不寻常的例子，但我们将其包含在内是为了举例说明具有适当先验的贝叶斯建模如何让我们修改和更改模型的几乎任何部分，前提是我们知道如何将我们的假设表示为模型过程的一部分。然后，我们使用蒙特卡洛方法转动推理曲柄并获得可靠的参数估计。引入更多符号，让 \(\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)\)，并在 \(\mathbf{s}_i\) 处而不是在\(\tilde{\mathbf{s}}_i\)。我们在这里选择的坐标系有些随意；我们还可以选择使用极坐标，并将先验放在位置误差的大小和角度上。在此示例中，比例 \(\sigma_s\) 被视为已知，以便模型表示不同的假设水平的坐标误差。 \[\begin{对齐} \Delta_i &\sim \mathrm{正态}(\mathbf{0}, \sigma_s^2\mathbf{I}_2) \\ \mathbf{s}_i &= \tilde{\mathbf{s}}_i + \Delta_i \\ \mu &\sim \mathrm{正态}(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 color from matplotlib.patches import Circle from stsp import use_sepia_gunmetal RANDOM_SEED = 8927 np 。随机的 。种子（RANDOM_SEED）主题，sepia_cmap = use_sepia_gunmetal（）plot_bg_color =“#1c1c1d”plot_text_color =“#e8e8e8”plt。风格 。使用 ("dark_background") mpl 。 rcParams 。更新（{“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",})主题["paper"]=plot_bg_color df=pd. read_csv（路径（“../../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 " ]