蒙特卡罗方法中的集合解与跨维度模型

# 蒙特卡罗方法中的集合解与跨维度模型 ## 1. 贝叶斯公式与后验概率 在相关计算中，贝叶斯公式是一个核心工具，其表达式为： \[p(m|d_{obs}) = \frac{p(d_{obs}|m)p_A(m)}{p(d_{obs})}\] 其中，\(p_A(m)\) 是先验概率，表达了对参数 \(m\) 值的先验期望；\(p(d_{obs}|m)\) 是似然概率，体现了在给定特定模型值时数据出现的概率；\(p(m) \equiv p(m|d_{obs})\) 是后验概率，也就是我们要求解的内容。在 Metropolis - Hastings 公式中，\(p(d_{obs})\) 这一因子会被消去，所以通常不再对其进行深入讨论。 先验概率和似然概率常常被假定为正态概率密度函数（pdf），具体形式如下： - 似然概率： \[p(d_{obs}|m) = (2\pi)^{-N/2} \det[\text{cov} d]^{-1/2} \exp[-E(m)]\] 其中 \(E(m) = e^T [\text{cov} d]^{-1} e\)，\(e = d_{obs} - g(m)\)。 - 先验概率： \[p_A(m) = (2\pi)^{-M/2} \det[\text{cov} m]_A^{-1/2} \exp[-L(m)]\] 其中 \(L(m) = l^T [\text{cov} m]_A^{-1} l\)，\(l = \langle m \rangle - m\)。 这里，\(E\) 是加权预测误差（相对于预测 \(g(m)\) 进行测量），\(L\) 是先验信息中的加权误差（相对于先验模型 \(\langle m \rangle\) 进行测量）。需要注意的是，\((2\pi)^{-N/2} \det[\text{cov} d]^{-1/2}\) 这一因子在 Metropolis - Hastings 公式中会被消去，因为数据的数量及其先验协方差都不会随 \(m\) 发生变化。 ## 2. Metropolis - Hastings 算法与后验概率的关系 Metropolis - Hastings 算法的作用仅仅是对后验概率密度函数 \(p(m|d_{obs})\) 进行采样，而并非创造它。后验概率是由使用者通过选择似然概率和先验概率密度函数来确定的（如上述公式所示）。一旦这些概率密度函数确定下来，后验概率的性质以及从其中抽取的集合的性质也就随之固定（存在一定的随机变化）。因此，在选择这些概率密度函数时需要格外谨慎，要确保它们能够体现数据的噪声特性（对于似然概率而言）以及先验信息（对于先验概率而言）。然而，在处理具有非线性模型和/或非正态概率密度函数的复杂问题时，这些选择所带来的后果可能难以预先估计，往往需要通过检查集合的性质才能发现。 ## 3. 集合解的应用 集合解（例如使用马尔可夫链蒙特卡罗（MCMC）方法生成的集合解）的用途不仅仅局限于推断单个估计解 \(m_{est}\)。实际上，我们可以从集合解中推断出 \(m_{est}\) 及其后验协方差 \([\text{cov} m]\)，具体方法是计算集合的样本均值和样本协方差。 一种可能的应用是利用 Backus - Gilbert 方法来探索模型分辨率和协方差之间的权衡关系。具体操作步骤如下： 1. 将问题视为 \(Gm = d_{obs}\) 形式的反问题，其中 \(G = I\)，\(d_{obs} = m_{est}\)，\([\text{cov} d] = [\text{cov} m]\)。 2. 该方程定义了与模型分辨率矩阵 \(R(\alpha = 1) = I\) 相关的唯一平均值，在权衡曲线上 \(\alpha = 1\) 这一点的协方差为 \(C_m(\alpha = 1) \equiv [\text{cov} m]\)。这里，\(\alpha\) 控制着模型分辨率的扩展和协方差大小的相对权重。 3. 利用 Backus - Gilbert 方法尝试找出一个不同的 \(\alpha\) 值，使得对应的局部平均值具有较差（但仍可接受）的模型分辨率和更好的方差。 ## 4. 集合解的示例 ### 4.1 示例 1：简单非线性反问题 我们对之前研究过的简单非线性反问题计算了一个包含一百万个成员的集合解。具体步骤如下： 1. 合成数据的生成：通过向真实数据添加均值为零、方差为 \(\sigma_d^2 = (0.3)^2\) 的随机噪声来创建合成数据。 2. 先验信息的添加：与之前的示例不同，我们加入了关于解接近先验值 \(m_A = [1.43, 1.50]^T\) 的先验信息，先验方差为 \(\sigma_m^2 = (0.5)^2\)。 3. 状态转移的计算：使用正态概率密度函数来计算 \(m\) 和 \(m'\) 之间的转移，公式为： \[q(m_i|m_i') = q(m_i'|m_i) = \frac{1}{\sqrt{2\pi} \sigma_q} \exp\left(-\frac{(m_i - m_i')^2}{2\sigma_q^2}\right)\] 我们使用 \(\sigma_q^2 = (0.05)^2\)，这样当前解和提议解之间的差异约为解可能范围的几个百分点。 4. 集合的计算：使用 Metropolis - Hastings 算法计算出一个包含一百万个成员的集合。这些成员在 \((m_1, m_2)\) 平面上形成一个紧密的簇，位于广义误差 \(\Psi = E + L\) 最小值的最深部分，其中心接近之前方法估计的点。 5. 预测数据的处理：计算每个成员模型的预测数据 \(d_i^{obs}\)，并进行分箱处理，以生成经验概率密度函数 \(p(d_i^{obs})\)。这些预测模型的离散程度比观测数据窄，并且在数据的峰值附近最窄，这表明数据中最受约束的属性是峰值位置。 ### 4.2 示例 2：拉普拉斯变换类问题 我们对一个拉普拉斯变换类问题计算了集合解，其中模型参数 \(m_j\) 与数据 \(d_i\) 之间的关系为： \[d_i = \sum_{j = 1}^{M} m_j \exp(-c_i z_j)\] 其中 \(z_j = \frac{j - 1}{M - 1}\)，\(M = 11\)，\(c\) 是范围在 0 到 0.6 之间的已知常数。真实模型是正弦波 \(m_{true}(z) = \sin(2\pi z)\)。具体步骤如下： 1. 合成数据的生成：通过向真实数据 \(d_{true}\)（通过将 \(m = m_{true}\) 代入上述公式得到）添加方差为 \(\sigma_d^2\) 的随机噪声来创建合成数据。 2. 先验概率的选择：先验概率 \(p_A(m)\) 被选择为正态、不相关、均值为零且方差为 \(\sigma_m^2 = (5)^2\) 的概率密度函数，这对应于较弱的先验信息，即解最有可能在 \(\pm 5\) 的范围内。 3. 状态转移概率的选择：\(q(m|m')\) 和 \(q(m'|m)\) 被选择为正态、不相关且方差为 \(\sigma_q^2\) 的概率密度函数（如示例 1 中的公式）。我们同样使用 \(\sigma_q^2 = (0.05)^2\)。 4. 集合的分析：集合的直方图显示，只有解的较浅（左侧）部分能够被准确恢复。这是拉普拉斯变换类问题的典型特征，即数据对模型的较深部分敏感度较低，因此较深部分主要受到先验信息的强烈影响（尽管先验信息较弱）。 5. 概率的评估：我们可以评估解具有某些有趣属性的概率。例如，通过确定集合中 \(m_2\) 是最大模型参数的模型的百分比（在这种情况下为 63%），可以评估模型在 \(z = 0.2\) 附近有峰值的概率。 ## 5. 跨维度模型的概念 想象你在公园散步时，看到空中有个东西呼啸而过，但你不确定它是什么。它可能是一个球，也可能是一只鸟。如果是球，那么用一个模型参数（例如它的直径 \(m_1\)）就足以描述它；但如果是鸟，则需要几个参数（例如它的长度 \(m_{10}\) 和翼展 \(m_{20}\)）。 这就是跨维度模型的世界！整体模型是由 \(D\) 种不同的模型类