Canoe Basis：Markov Chain Mont Carlo

It has been 501 days since the last update, the content of the article may be outdated.

MCMC方法

贝叶斯统计学基础

在贝叶斯统计学中，我们试图估计未知参数（用 θ 表示），给定观测到的数据（用 D 表示）。这个过程涉及到两个主要的概率分布：

1. 先验分布（Prior Distribution）：表示在观测到任何数据之前我们对参数的信仰。通常用 P(θ) 表示。

2. 似然函数（Likelihood Function）：表示给定参数情况下观测到数据的概率。通常用 P(D|θ) 表示。

这两者结合起来构成了贝叶斯定理：

$$P(\theta |D)=\frac{P(D|\theta )\cdot P(\theta )}{P(D)}$$

其中，

• $P(\theta |D)$ 是后验概率，即在观测到数据后参数的概率分布；
• $P(D|\theta )$ 是似然函数；
• $P(\theta )$ 是先验概率；
• $P(D)$ 是数据的边际概率，也就是在所有可能的参数下观测到这个数据的概率。

---

马尔可夫链

马尔可夫链 是一种随机过程，根据一定的概率规则在一系列状态之间进行转移。马尔可夫链的定义特征是马尔可夫性质，即未来状态只依赖于当前状态，而不依赖于先前事件的序列。这是一个重要的性质，因为MCMC利用了这个性质来构建一个接受-拒绝的采样算法。
马尔可夫链广泛应用于物理学、经济学、计算机科学和生物学等领域，用于建模表现出无记忆性的随机过程和系统。

矩阵表示

马尔可夫链的演变可以使用转移概率矩阵进行简洁的表示。假设我们的马尔可夫链有 $N$ 个状态，从状态 $i$ 转移到状态 $j$ 的概率记为 $P_{ij}$。转移概率矩阵 $P$ 则表示为：
$$\begin{array}{r}P=\left[\begin{array}{cccc}{P}_{11}& P_{12}& \cdots & P_{1N}\\ P_{21}& P_{22}& \cdots & P_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ P_{N1}& P_{N2}& \cdots & P_{NN}\end{array}\right]\end{array}$$
这里，每一行代表当前状态，每一列代表下一个状态。元素 $P_{ij}$ 表示从状态 $i$ 转移到状态 $j$ 的概率。

[图片来源:【Nature Methods】Markov models—Markov chains](https://www.nature.com/articles/s41592-019-0476-x)

精细平衡和守恒条件

为了使马尔可夫链达到稳态分布（随时间不变的分布），通常需要满足**详细平衡(detailed balance)**条件。

• 精细平衡：

若对于所有状态 $i$ 和 $j$，以下方程成立，则马尔可夫链满足详细平衡：

$${\pi }_i\cdot P_{ij}={\pi }_j\cdot P_{ji}$$

这里，${\pi }_i$ 表示处于状态 $i$ 的稳态概率。详细平衡确保链处于平衡状态，概率不会在一个方向上累积。

• 守恒条件：

1. 不可约性：若在有限步内可以从任一状态到达任何其他状态，则马尔可夫链是不可约的。这确保没有孤立的状态集。

2. 非周期性：若转移没有固定的模式，则马尔可夫链是非周期的。这意味着链不会陷入循环，所有状态都有可能在任何时刻被重新访问。

满足这些条件确保了马尔可夫链的长期行为收敛到唯一的稳态分布。

总结而言，马尔可夫链是一个数学模型，描述系统通过一系列状态的转移概率。转移过程通常用转移概率矩阵表示，详细平衡条件确保存在一个稳态分布。不可约性和非周期性是链收敛到稳态分布的重要性质。

---

Metropolis-Hastings 算法

Metropolis-Hastings 算法是一种基于马尔可夫链的 MCMC 抽样方法，它通过接受-拒绝的机制在参数空间中游走，从而生成参数的后验分布的样本。

算法步骤：

1. 初始化：选择初始参数值 ${\theta }_0$。

2. 提议分布（Proposal Distribution）：选择一个用于提议新参数值的分布，通常表示为 $q({\theta }^{\prime }|\theta )$。

3. 迭代：对于每个迭代 $t$，重复以下步骤：

   a. 从提议分布中抽取一个候选参数值 ${\theta }^{\prime }$。

   b. 计算接受比率（Acceptance Ratio）：
   $$r=\frac{P({\theta }^{\prime }|D)}{P({\theta }_t|D)}\cdot \frac{q({\theta }_t|{\theta }^{\prime })}{q({\theta }^{\prime }|{\theta }_t)}$$

   c. 以概率 $min(r,1)$ 接受候选参数值：
   $${\theta }_{t+1}=\{\begin{array}{ll}{\theta }^{\prime }& \text{with probability }min(r,1)\\ {\theta }_t& \text{with probability }1-min(r,1)\end{array}$$

4. 重复步骤 3 直到获得足够的样本。

统计学基础：

• 接受率（Acceptance Rate）：表示接受候选参数值的比率。合适的接受率通常在 0.2 到 0.5 之间。

• 遍历性（Ergodicity）：Metropolis-Hastings 算法需要保证马尔可夫链是遍历的，即在足够的迭代次数后，任何状态都有可能被访问。

• 平稳性（Stationarity）：Metropolis-Hastings 算法在达到平稳状态后生成的样本服从目标后验分布。

• 收敛性（Convergence）：检查算法是否收敛到目标分布。一种常见的方法是通过观察连续的样本来检查参数估计值是否趋于稳定。

• Burn-in 阶段：初始迭代可能由于初始参数的选择而不够准确，因此可能需要丢弃一些初始样本（称为 Burn-in 阶段）。

请注意，以上的方程和描述是对 MCMC 的高层次概括，实际应用中需要根据具体问题进行调整和细化。

---

MCMC采样器

下面是一个使用Metropolis-Hastings MCMC方案实现符合给定概率分布样本抽样的C++一般结构：

#include <iostream>
#include <vector>
#include <cstdlib>   // For rand() and RAND_MAX
#include <cmath>     // For exp()

// Abstract likelihood function
double likelihood(const std::vector<double>& data, const std::vector<double>& parameters) {
    // Calculate the likelihood of the data given the parameters
    // ...
}

// Abstract prior function
double prior(const std::vector<double>& parameters) {
    // Calculate the prior probability of the parameters
    // ...
}

// Abstract posterior function
double posterior(const std::vector<double>& data, const std::vector<double>& parameters) {
    // Calculate the posterior probability of the parameters
    return likelihood(data, parameters) * prior(parameters);
}

// Abstract proposal function
std::vector<double> proposeNewParameters(const std::vector<double>& currentParameters) {
    // Propose new parameters from a proposal distribution
    // ...
}

// Abstract Metropolis-Hastings MCMC function
std::vector<std::vector<double>> metropolisHastings(const std::vector<double>& data,
                                                    const std::vector<double>& initialParameters,
                                                    int numIterations) {
    std::vector<double> currentParameters = initialParameters;
    std::vector<std::vector<double>> samples;

    for (int i = 0; i < numIterations; ++i) {
        // Propose new parameters from a proposal distribution
        std::vector<double> proposedParameters = proposeNewParameters(currentParameters);

        // Compute acceptance ratio
        double acceptanceRatio = std::min(1.0, posterior(data, proposedParameters) / posterior(data, currentParameters));

        // Accept or reject the proposed parameters
        if (static_cast<double>(rand()) / RAND_MAX < acceptanceRatio) {
            currentParameters = proposedParameters;
        }

        // Record samples
        samples.push_back(currentParameters);
    }

    return samples;
}

int main() {
    // Example usage
    std::vector<double> data = { /* ... your observed data ... */ };
    std::vector<double> initialParameters = { /* ... initial guess for parameters ... */ };
    int numIterations = /* ... number of MCMC iterations ... */;

    // Run MCMC
    std::vector<std::vector<double>> samples = metropolisHastings(data, initialParameters, numIterations);

    // Analyze and interpret the samples to obtain parameter estimates
    // ...

    return 0;
}

---

MCMC参数反演

问题描述：现有大气层顶（TOA）的微波亮温观测 $T{B}_{obs}$ 和辐射传输模型，可以使用{Temperature，RelativeHumidity}等参数来模拟TOA亮温 $T{B}_{simu}$ 。请使用MCMC方法，根据现有观测和模式，进行大气双参数的反演

从大气层顶部 (TOA) 亮度温度观测中检索温度和相对湿度等大气参数涉及一个复杂的过程，通常需要马尔可夫链蒙特卡罗 (MCMC) 等复杂方法。下面是使用通用抽象语言的简化示例。请注意，实际实现可能会根据您的问题的具体细节而有所不同，并且您可能需要根据您的需求调整代码。

#include <iostream>
#include <vector>
#include <random>

// Function to simulate TOA brightness temperature based on atmospheric parameters
double simulateTOA(const double temperature, const double relativeHumidity) {
    // Your forward model implementation goes here
    // This is a simplified example; replace it with your actual forward model
    double simulatedTB = temperature + relativeHumidity + 0.1 * rand() / RAND_MAX;
    return simulatedTB;
}

// Function to calculate the likelihood of observed TOA brightness temperature given simulated TOA brightness temperature
double likelihood(const double obsTB, const double simuTB, const double sigma) {
    // Assuming a simple Gaussian likelihood function
    return exp(-(obsTB - simuTB) * (obsTB - simuTB) / (2 * sigma * sigma));
}

// Function to perform Metropolis-Hastings MCMC
void MCMC(const std::vector<double>& obsTB, double& estimatedTemperature, double& estimatedRH) {
    const int iterations = 10000;
    const double proposalSigma = 0.1;

    // Initialize parameters
    double currentTemperature = 288.0; // Initial guess for temperature
    double currentRH = 0.5; // Initial guess for relative humidity

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> uniformDist(0.0, 1.0);

    for (int i = 0; i < iterations; ++i) {
        // Propose new parameters
        double proposedTemperature = currentTemperature + proposalSigma * (2 * uniformDist(gen) - 1);
        double proposedRH = currentRH + proposalSigma * (2 * uniformDist(gen) - 1);

        // Simulate TOA brightness temperature with proposed parameters
        double simuTB = simulateTOA(proposedTemperature, proposedRH);

        // Calculate likelihoods
        double likelihoodCurrent = 1.0;
        double likelihoodProposed = 1.0;
        for (size_t j = 0; j < obsTB.size(); ++j) {
            likelihoodCurrent *= likelihood(obsTB[j], simulateTOA(currentTemperature, currentRH), 0.1);
            likelihoodProposed *= likelihood(obsTB[j], simuTB, 0.1);
        }

        // Accept or reject the proposed parameters
        double acceptanceRatio = likelihoodProposed / likelihoodCurrent;
        if (uniformDist(gen) < acceptanceRatio) {
            currentTemperature = proposedTemperature;
            currentRH = proposedRH;
        }
    }

    // Set the estimated values
    estimatedTemperature = currentTemperature;
    estimatedRH = currentRH;
}

int main() {
    // Example usage
    std::vector<double> observedTOA = {288.1, 288.2, 287.8}; // Replace with your actual observed TOA brightness temperatures

    double estimatedTemperature, estimatedRH;
    MCMC(observedTOA, estimatedTemperature, estimatedRH);

    std::cout << "Estimated Temperature: " << estimatedTemperature << std::endl;
    std::cout << "Estimated Relative Humidity: " << estimatedRH << std::endl;

    return 0;
}

此示例使用带有高斯似然函数的简单 Metropolis-Hastings MCMC 算法。将函数替换simulateTOA为您的实际正向模型。此外，根据问题的特征调整似然函数和其他参数。

问题

1. 如何形象地理解MCMC方法walker过程探索概率分布？

   为了形象地理解马尔可夫链蒙特卡洛（MCMC）方法中的”walk”过程是如何构建概率分布的，我们可以想象在参数空间中漫步的过程。以下是一个简单的比喻：

   比喻
   想象你身处一个大草原，你的目标是找到一个藏宝点（概率分布的峰值），而你对草原的地形一无所知。你手里有一张地图，但是地图只显示你当前所在位置的周围区域。

   步骤:1.初始位置

   • 你从草原的某个位置开始，这是你的初始猜测。
     2. 漫步过程

   • 你开始漫步，每一步都是在当前位置的基础上随机选择一个方向走。

   • 每一步的大小和方向都是根据一个特定的规则（转移概率分布）确定的，这就是马尔可夫链的特性。
     3. 接受或拒绝

   • 到达一个新位置后，你会考虑这个新位置相比于当前位置的好坏（通过目标概率分布来评估）。
     如果新位置更有利于找到藏宝点，你会接受这个新位置，否则你可能以一定概率拒绝它。

   • 这一步保证了在参数空间中朝着更有可能找到藏宝点的方向前进。
     4. 反复迭代

   • 你反复执行上述步骤，每一步都根据当前位置和转移概率分布进行更新。

   • 随着漫步的进行，你的位置逐渐集中在藏宝点附近，这就是MCMC方法构建概率分布的过程。

   关键点:

   1. 漫步的起始点：初始位置通常是在参数空间中随机选择的，代表了对概率分布的初始猜测。

   2. 转移概率分布：决定了每一步漫步的方向和大小，是马尔可夫链的关键。

   3. 目标概率分布：用于评估新位置的好坏，通常是我们希望找到的概率分布。

   4. 接受或拒绝规则：确保朝着概率分布的峰值前进，而不是陷入局部最小值。

   通过这个比喻，你可以想象在参数空间中通过随机漫步的过程中，逐渐找到概率分布的峰值，这就是MCMC方法的核心思想。这种随机漫步的方式使得我们能够探索整个参数空间，而不仅仅是依赖于确定性的搜索方法。

2. 新位置不好的情况下，我为什么还要有一定概率接受它呢？这是为了满足什么条件

   接受新位置的概率是为了确保马尔可夫链蒙特卡洛（MCMC）方法在参数空间中能够全面地探索，而不会陷入局部最小值或局部最大值。这涉及到MCMC方法中的一个重要概念，即Metropolis-Hastings算法，其中包含了接受或拒绝的机制，以满足精细平衡条件。

   在MCMC方法中，我们希望构建一个马尔可夫链，使其平稳分布（稳态分布）与我们想要采样的目标概率分布一致。这就是精细平衡条件，具体表达为：
   $$\begin{array}{}\text{(1)}& \pi (x)\cdot T(x\to x^{\prime })=\pi (x^{\prime })\cdot T(x^{\prime }\to x)\end{array}$$
   其中，

   • $\pi (x)$ 是我们想要采样的目标概率分布。
   • $T(x\to x^{\prime })$ 为从 $x$ 状态到 $x^{\prime }$ 状态地转移概率。

   Metropolis-Hastings算法是一种基于精细平衡条件的MCMC抽样方法，接受或拒绝：以概率min(r,1) 接受新状态，否则维持在当前状态。这保证了在新状态较好时接受，而在新状态较差时以一定概率也接受，以便在整个参数空间中进行探索。
   通过这个机制，即使新状态相对较差，也有一定的概率接受它。这是为了确保全局探索，从而不会被困在局部最小值或局部最大值。在长时间运行中，Metropolis-Hastings算法使得马尔可夫链趋向于平稳分布，即目标概率分布。

3. 如何理解精细平衡？

   精细平衡条件确保了在平稳状态下，系统中的每一对状态之间的转移速率（或概率）是相等的。换句话说，系统在这些状态之间保持平衡，不会有概率流向一个方向而不是另一个方向。

4. 电子能级跃迁的精细平衡

   考虑一个原子系统，其中有两个能级：基态 $E_1$ 和激发态 $E_2$ 。电子可以从基态跃迁到激发态，也可以从激发态返回基态。这些跃迁过程可以通过发射或吸收光子来实现。

   精细平衡条件
   精细平衡条件要求在平稳状态下，电子从基态到激发态的跃迁速率与电子从激发态返回基态的跃迁速率之间存在平衡。具体表达为：
   $$N_1\cdot A_{12}=N_2\cdot A_{21}$$

   其中，

   • $N_1$ 和 $N_2$ 分别是基态和激发态的电子数。
   • $A_{12}$ 是从基态到激发态的跃迁速率。
   • $A_{21}$ 是从激发态返回基态的跃迁速率。
     在这个过程中，光子的吸收和发射实际上是驱动电子能级之间跃迁的过程。如果精细平衡条件得到满足，系统将达到一个稳定的状态，其中电子数在不同能级之间保持平衡。
     也就是说，就算激发态的电子很少（激发概率小），但也是有一定概率跃迁到激发状态的，如果这部分概率被忽略了，两个状态的电子数目就无法达到平衡。

5. Why and how to use logprob in MCMC reference?

   数值稳定性：

   1. 避免下溢：

      • 概率可以非常小，特别是在处理连续数据或高维参数空间时。
      • 许多小概率相乘可能导致数值下溢，即值变得太小以至于计算机无法准确表示。

   2. 对数变换：

      • 取概率的对数将乘法操作转换为加法操作。
      • 对数变换有助于防止数值下溢，使计算更稳定。

   简化计算：

   3. 加法代替乘法：

      • 在MCMC中，计算通常涉及乘以似然和先验概率。
      • 取对数将这些乘积转换为和，简化了计算。

   4. 在对数空间中的求和：

      • 在计算后验概率时，以对数空间工作允许有效地进行求和而不是乘法。

   MCMC算法：

   5. 提案接受概率：
      • MCMC算法涉及提出新的参数值，并根据接受概率接受或拒绝这些提议。
      • 这个接受概率通常用后验概率的比率来表示，其中涉及对数概率。

   示例伪代码：

   # 计算对数似然和对数先验
   log_likelihood = calculate_log_likelihood(parameters, data)
   log_prior = calculate_log_prior(parameters)
   
   # 合并成对数后验
   log_posterior = log_likelihood + log_prior
   
   # 在MCMC提议步骤中，使用对数概率计算接受比率
   acceptance_ratio = min(1, exp(log_posterior - current_log_posterior))
   
   # 如果提议被接受，则用提议状态更新当前状态
   if random_uniform(0, 1) < acceptance_ratio:
       current_parameters = proposed_parameters
       current_log_posterior = log_posterior
   
       

最后附上成哥给我上课的白板：

MCMC讲解（局部）