[Machine Learning] 采样

---

蒙特卡罗方法 (Monte Carlo Methods)

逆变换法 (Inverse Transform Sampling)

逆变换法的基本思想是，如果U是[0,1]上的均匀分布，那么累积分布函数（CDF）的逆函数应用到U上将产生符合目标分布的随机样本。

假设我们有一个随机变量X，其概率密度函数（PDF）为f(x)，累积分布函数为F(x)。那么，F(x)的逆函数F^(-1)(u)，应用于均匀分布的随机变量U，将产生符合目标分布f(x)的随机样本。

具体步骤如下：

• 生成一个[0,1]上的均匀分布的随机数u。
• 计算F^(-1)(u)，得到的结果就是一个符合目标分布f(x)的随机样本。

注意，F^(-1)(u)并不总是容易计算，因此这种方法并不总是可用。

拒绝采样法 (Rejection Sampling)

拒绝采样的原理是通过生成一个已知分布（包络线分布）的样本集，然后基于目标分布和包络线分布的概率密度比例来决定是否接受这些样本点，最终得到一个近似于目标分布的样本集。这个方法的基本原理可以从接受-拒绝准则（Accept-Reject Criterion）来理解。

接受-拒绝准则的基本思想是，如果我们能够找到一个概率密度函数为g(x)的分布（即包络线分布），使得对于所有的x，都有c * g(x) >= f(x)，其中f(x)是目标分布的概率密度函数，c是一个大于等于1的常数。那么我们可以通过以下步骤生成目标分布的样本：

1. 从已知分布g(x)中采样一个点x。
2. 从均匀分布U(0,1)中采样一个随机数u。
3. 计算接受概率：r = f(x) / (c * g(x))。
4. 如果u <= r，则接受这个样本点x；否则，拒绝这个样本点，并重复步骤1。

这个过程背后的直观原理是，我们尝试用一个简单的已知分布g(x)来近似目标分布f(x)。在这个过程中，我们生成了很多样本点，但是并不是所有的样本点都能被接受。实际上，只有当一个样本点x的概率密度在目标分布f(x)中相对较高时，它才会被接受。换句话说，接受概率r度量了样本点x在目标分布f(x)中的“重要性”。通过这种方式，拒绝采样方法可以生成一个与目标分布f(x)相近的样本集。

需要注意的是，拒绝采样方法的效率很大程度上取决于包络线分布g(x)与目标分布f(x)的匹配程度。如果g(x)能够很好地近似f(x)，那么拒绝率会较低，采样效率较高。相反，如果g(x)与f(x)匹配得不好，那么拒绝率会较高，采样效率较低。因此，在实际应用中，选择一个合适的包络线分布g(x)是很重要的。

重要性采样法 (Importance Sampling)

假设我们有一个难以直接采样的目标分布$p(x)$，我们想要估计的是某个函数$f(x)$关于$p(x)$的期望，即：

$$E_p[f]=\int f(x)p(x)dx$$

由于我们无法直接从$p(x)$中采样，我们选择一个易于采样的重要性分布$q(x)$，然后我们可以将上面的期望写成：

$$E_p[f]=\int f(x)\ast [p(x)/q(x)]\ast q(x)dx=E_q[f\ast (p/q)]$$

其中，$E_q\left[f\frac{p}{q}\right]$表示的是函数$f(x)\frac{p(x)}{q(x)}$关于$q(x)$的期望。

这样，我们就可以通过以下步骤进行重要性采样：

1. 从重要性分布$q(x)$中抽取$N$个样本${x^{(i)}}_{i=1}^N$。
2. 计算每个样本的权重$w^{(i)}=\frac{p(x^{(i)})}{q(x^{(i)})}$，这些权重用于校正从$q(x)$采样带来的偏差。
3. 根据权重和函数值计算加权平均作为期望的估计值，即：$\hat{E}p[f]=\frac{1}{N}\sum {i=1}^N{w}^{(i)}f(x^{(i)})$。

---

马尔可夫链蒙特卡罗方法 (Markov Chain Monte Carlo, MCMC)

For Bayesian statistics, the model parameters has a prior.

$$p(\theta |D)=\frac{p(D|\theta )\cdot p(\theta )}{\int p(D|{\theta }^{\prime })\cdot p({\theta }^{\prime })d{\theta }^{\prime }}$$

• Generally, it’s very difficult to exactly compute $p(\theta |D)$
• Sample from $p(\theta |D)$: using MCMC
• Idea: construct a Markov chain such that its stationary distribution is the desired density function

马尔可夫链 (Markov Chains)

• Markov chain is a stochastic process that follows the Markov property.
• Markov property means that the future state of the process only depends on the current state.
• Almost memory-less system.

转移概率矩阵 (Transition Probability Matrix)

对于具有 $n$ 个状态的马尔可夫链，其转移概率矩阵 $P$ 是一个 $n\times n$ 的矩阵，其中每个元素 $p_{ij}$ 表示在任意时刻，从状态 $i$ 转移到状态 $j$ 的概率，即 $p_{ij}=P(X_{t+1}=j|X_t=i)$。由于每个状态都必须转移到某个状态（包括可能的自我转移），所以转移概率矩阵的每一行之和都等于1。

不变分布 (Invariant Distribution) / 平稳分布 (Stationary Distribution)

在马尔可夫链中，如果存在一个概率分布$\pi$，使得对于转移概率矩阵$P$，有$\pi P=\pi$，那么我们就称$\pi$为该马尔可夫链的不变分布或平稳分布。

更准确地说，如果我们有一个离散状态空间的马尔可夫链，其状态空间为$S$，转移概率矩阵为$P$，那么一个概率分布$\pi :S\to [0,1]$被称为该马尔可夫链的不变分布或平稳分布，如果它满足以下条件：

$${\pi }_i=\sum _{j\in S}{\pi }_j{P}_{ji}\forall i\in S$$

这个条件表示，如果马尔可夫链在时间$t$的状态分布为$\pi$，那么在时间$t+1$的状态分布仍然是$\pi$。

在这里，$P_{ji}$是从状态$j$转移到状态$i$的概率，${\pi }_i$是马尔可夫链在状态$i$处的概率，$S$是所有可能的状态。这个公式可以理解为，无论马尔可夫链在任何一个状态，其下一个状态的分布都是不变的，这就是为什么我们称之为"不变"或"平稳"分布的原因。

Eventually converges to an invariant distribution

• 马尔可夫链必须是不可约的（irreducible），即不存在状态之间的隔离。这意味着任意两个状态之间存在一条路径，使得在有限步数内可以从一个状态转移到另一个状态。

• 马尔可夫链必须是非周期的（aperiodic），即不存在周期性的状态循环。这确保了系统没有固定的时间模式，从而使得长期行为更加随机和无序。

或者说，马尔可夫链是遍历的（ergodic），即从任意状态出发，经过一定的步数后可以到达任何其他状态。这保证了系统在有限时间内可以从任意状态转移到任意其他状态。

在满足这些条件的前提下，马尔可夫链的不变分布可以通过迭代转移矩阵的操作得到。当马尔可夫链的状态转移达到稳定状态时，即使经过多次迭代，状态分布也不再发生显著变化，这就是马尔可夫链收敛到不变分布的过程。

这表示马尔可夫链的长期行为是可预测的，并且与初始状态无关。

证明

假设我们有一个马尔可夫链，其转移概率矩阵 $P$ 是一个随机矩阵，且所有的状态都是连通的并且是非周期性的。根据 Perron-Frobenius 定理，我们知道这样的矩阵存在一个且只有一个主特征值，且主特征值的绝对值等于 1，并且对应的特征向量是正的。我们将主特征值记为 ${\lambda }_1$，对应的特征向量记为 $v_1$，满足：

$$P{v}_1={\lambda }_1{v}_1$$

由于 ${\lambda }_1$ 的绝对值等于 1，我们可以设 ${\lambda }_1=1$。此时，我们有：

$$P{v}_1=v_1$$

这表明 $v_1$ 是一个平稳分布。

至于马尔可夫链最终会收敛到平稳分布，我们可以通过对转移概率矩阵 $P$ 的连续迭代来理解。对于任意的初始分布 $p_0$，我们有：

$$p_{t+1}=p_{t}P$$

通过连续迭代，我们有：

$$p_t=p_0{P}^t$$

当 $t$ 趋向于无穷大时，由于 $P$ 的其他特征值的绝对值都小于 1，因此 $P^t$ 的影响将主要由主特征值决定，即 $P^t$ 将趋向于一个矩阵，其每一行都等于 $v_1$。因此，$p_t$ 将收敛到平稳分布 $v_1$。

细致平衡条件 (Detailed Balance Condition)

细致平衡条件是马尔可夫链达到平稳分布的一个充分条件，但并不是必要条件。对于一个马尔可夫链，如果存在一个概率分布 $\pi$，使得对所有的状态 $i$ 和 $j$，以下等式成立：

$$\pi (i)\cdot p_{ij}=\pi (j)\cdot p_{ji}$$
那么，我们就说这个马尔可夫链满足细致平衡条件。

简单来说，细致平衡条件就是在平稳状态下，从状态 $i$ 转移到状态 $j$ 的总概率等于从状态 $j$ 转移到状态 $i$ 的总概率。

Metropolis-Hastings 算法

MH算法是一种常用的MCMC方法，它的目标是生成一个样本序列，这个序列的分布可以逼近我们感兴趣的目标分布。基本步骤如下：

1. 选择一个初始状态 $z^{(0)}$ 和一个提议分布 $Q(z^{\prime }|z)$，其中 $Q(z^{\prime }|z)$ 表示在给定当前状态 $z$ 的情况下生成新状态 $z^{\prime }$ 的概率。
   
   

2. 对于 $i=0,1,2,\dots$ 执行以下步骤：

• 从提议分布 $Q(z^{\prime }|z^{(i)})$ 中抽取一个候选状态 $z^{\prime }$。

• 计算接受概率 $\alpha =\min \left\{1,\frac{p(z^{\prime })Q(z^{(i)}|z^{\prime })}{p(z^{(i)})Q(z^{\prime }|z^{(i)})}\right\}$。
  $p(z)$ 通常是我们想要进行采样的目标分布，通常我们会有它的表达式或者至少能计算其概率密度函数在某个点的值。$p(z^{\prime })$ 和 $p(z^{(i)})$ 分别就是目标分布在新状态 $z^{\prime }$ 和当前状态 $z^{(i)}$ 处的概率密度。 $Q(z^{(i)}|z^{\prime })$ 和 $Q(z^{\prime }|z^{(i)})$ 则是提议分布在相应点的概率密度。

• 从均匀分布 $U(0,1)$ 中抽取一个随机数 $u$。

• 如果 $u<\alpha$，则接受新的候选状态，即设置 $z^{(i+1)}=z^{\prime }$。否则，保持当前状态，即设置 $z^{(i+1)}=z^{(i)}$。

这样，通过反复迭代，我们就能生成一个样本序列，这个序列的分布会逐渐逼近目标分布 $p(z)$。

Symmetric Proposal Distribution

在某些情况下，我们可以选择一个对称的提议分布，这意味着从 $z$ 跳转到 $z^{\prime }$ 的概率与从 $z^{\prime }$ 跳转到 $z$ 的概率是相同的，即 $Q(z^{\prime }|z)=Q(z|z^{\prime })$。在这种情况下，提议分布是对称的。

在 Metropolis-Hastings 算法中，我们需要决定是否接受提议的状态转移。这是通过计算接受概率来实现的。当我们的提议分布是对称的时，接受概率的计算可以简化为：

$$\alpha =\min \left\{1,\frac{p(z^{\prime })}{p(z)}\right\}$$

这里，$p(z)$ 是我们想要采样的目标分布。

一个常见的对称提议分布的选择是正态分布，这意味着新的状态是在当前状态的基础上加上一个服从正态分布的随机量，即 $z^{\prime }\sim N(z,\sigma )$。

在这种情况下，当前的参数值 $\theta$ 是正态分布的均值，而标准差 $\sigma$ 则决定了新的参数值 ${\theta }^{\prime }$ 的范围。这就好像我们从当前状态 $\theta$ 转移到新的状态 ${\theta }^{\prime }$，并且转移的可能性由正态分布给出。

当我们使用MH算法时，我们的目标是找到一个状态转移矩阵（或者说提议分布），使得我们从任何初始状态开始，通过多次迭代，最终得到的状态的分布趋近于我们想要的目标分布 $p(\theta |D)$。这就是为什么我们需要计算接受率 $\alpha$，并决定是否接受新的参数值：我们希望保证算法最终收敛到正确的分布。

$$\alpha =\min \left\{1,\frac{p(z^{\prime })Q(z^{(i)}|z^{\prime })}{p(z^{(i)})Q(z^{\prime }|z^{(i)})}\right\}$$

而如果我们选择正态分布作为提议分布，那么这些值可以通过正态分布的概率密度函数计算得到。假设我们选择的正态分布以当前状态 $z^{(i)}$ 为均值，以 $\sigma$ 为标准差，那么 $Q(z^{\prime }|z^{(i)})$ 就是以 $z^{(i)}$ 为均值，$\sigma$ 为标准差的正态分布在 $z^{\prime }$ 处的概率密度。

$$Q(z^{\prime }|z^{(i)})=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(z^{\prime }-z^{(i)}{)}^2}{2{\sigma }^2}\right)$$

同样，$Q(z^{(i)}|z^{\prime })$ 就是以 $z^{\prime }$ 为均值，$\sigma$ 为标准差的正态分布在 $z^{(i)}$ 处的概率密度。

需要注意的是，这里的 $\sigma$ 会影响算法的表现。如果 $\sigma$ 太大，那么每一步的跳跃会很大，可能会导致接受率很低；如果 $\sigma$ 太小，那么每一步的跳跃会很小，可能需要很多步才能充分探索状态空间。通常，$\sigma$ 的选择需要根据具体问题来调整。

Example

对于抛硬币，观察了一系列的硬币投掷结果，并希望估计出硬币正面朝上的概率 $\theta$。这是一个参数估计的问题，但是我们需要对这个参数 $\theta$ 有一个先验的认识。在这个例子中，$\theta$ 是正面朝上的概率。

在贝叶斯学派的观点中，$\theta$ 是一个随机变量。我们可以计算给定观察到的数据后 $\theta$ 的后验分布。我们首先需要一个先验分布 $p(\theta )$，它代表了我们在看到数据之前对 $\theta$ 的认识。例如，我们可以假设 $\theta$ 服从均匀分布，即在 $0$ 到 $1$ 之间的所有值都是等可能的。然后，我们可以使用观察到的数据来计算似然函数 $p(D|\theta )$，这是硬币投掷结果的概率。

在这个例子中，我们可以使用 MCMC 算法来估计后验分布 $p(\theta |D)$。首先，我们选择一个初始点 ${\theta }_0$，然后迭代地生成新的 $\theta$ 值。在每一步，我们可以使用一个建议分布 $Q({\theta }^{\prime }|\theta )$ 来生成一个新的 ${\theta }^{\prime }$，然后使用 MH算法中的接受概率来决定是否接受这个新的值。

如果接受新的值，那么 ${\theta }^{\prime }$ 成为下一步的当前值；如果不接受，那么当前值保持不变。通过这个过程，我们可以生成一系列的样本，这些样本的分布近似于后验分布 $p(\theta |D)$。这样，我们就可以得到 $\theta$ 的一个分布，而不仅仅是一个点估计。

在这个例子中，$p(D|\theta )$ 可以由二项分布给出，因为每次投掷硬币都是独立的，且正面朝上的概率为 $\theta$。因此，给定 $\theta$，观察到 $k$ 次正面和 $n-k$ 次反面的概率为

$$p(D|\theta )=(\genfrac{}{}{0ex}{}{n}{k}){\theta }^k(1-\theta {)}^{n-k}$$

其中，$n$ 是投掷次数，$k$ 是正面朝上的次数，$\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ 是组合数，表示从 $n$ 次投掷中选择 $k$ 次正面的方式数。

对于先验分布 $p(\theta )$，我们可以选择一个简单的分布，比如 Beta 分布。Beta 分布在 $[0,1]$ 范围内定义，因此适合描述概率。对于参数为 1 的 Beta 分布，它是均匀分布，这表示我们在没有任何信息的情况下认为所有的 $\theta$ 都是等可能的。

一旦我们有了先验分布和似然函数，我们就可以使用 MCMC 算法来生成后验分布的样本。在每一步，我们可以使用一个建议分布，比如正态分布，来生成新的 ${\theta }^{\prime }$，然后计算接受率 $\alpha$，并决定是否接受这个新的值。通过多次迭代，我们可以得到一系列的 $\theta$ 样本，它们的分布近似于后验分布 $p(\theta |D)$。

在实践中，我们通常会在一开始的几步（被称为 “burn-in” 阶段）后，才开始从 MCMC 生成的样本中进行采样，以确保我们的样本来自平稳分布。我们还需要确保 MCMC 链的混合时间足够长，以得到独立的样本。通过分析这些样本，我们可以得到后验分布的各种特性，如均值、中位数、置信区间等。

MCMC practical considerations

• The samples at the start of the MCMC, before the algorithm converges to the true distribution are known as the burn-in
  period

  • It should be discarded
    
    

• The samples generated by MCMC are correlated since they are from a Markov chain

  • Previously, many practitioners advocated thinning the samples by taking say every $k^{th}$ sample
  • This was done for a few reasons historically
    • Reduce correlations and compute standard errors more easily
    • Less space needed to store the chain