1. 接受-拒绝重采样方法

重采样方法由两个步骤组成，第一个步骤提供近似分布的随机量，第二个步骤是校正机制。 在下文中，我们用 $f$ 表示目标分布，用 $g$ 表示辅助分布。

本课程我们只研究经典的接受-拒绝方法。

我们的目标是从概率密度为$f$的分布中抽取一个随机量$X$，这可能很难抽取。 我们引入辅助概率分布$g$并采取以下假设：
(A) 假设存在一些 $M>0$ 使得
$$\frac{f(x)}{g(x)}\le M\text{ 对于所有 }x\text{。 }$$

接受-拒绝方法的伪代码：

• 从分布 $g$ 中采样随机数量 Y
• 从 $U([0,1])$ 中抽取随机数量 $U$。
• • 如果 $U⩽\frac{f(Y)}{Mg(Y)}$，接受，即。 $X=Y$
• • 如果 $U>\frac{f(Y)}{Mg(Y)}$，则拒绝并重复该过程。

该算法的理论证明：

让我们考虑一维情况。 对于任何 $x\in \mathbb{R}$，考虑
$$\begin{array}{rl}P(X⩽x)& =P\left(Y⩽x\mid \frac{f(Y)}{Mg(Y)}⩾U\right)\\ & =\frac{P\left(Y⩽x,\frac{f(Y)}{Mg(Y)}⩾U\right)}{P\left(\frac{f(Y)}{Mg(Y)}⩽U\right)}\\ & =\frac{{\int }_{-\infty }^x\left({\int }_0^{\frac{f(Y)}{Mg(y)}}du\right)g(y)dy}{{\int }_{-\infty }^{+\infty }\left({\int }_0^{\frac{f(y)}{Mg(y)}}du\right)g(y)dy}\left(\because fracf(y)Mg(y)⩽1\text{ 对于所有 }y\right)\\ & =\frac{{\int }_{-\infty }^{x}f(y)dy}{{\int }_{-\infty }^{+\infty }f(y)dy}\end{array}$$
因此，$X$ 具有密度为 $f$ 的分布。

2. Python编程实现

接受-预测方法可用于对分布乘上一个常数的随机量进行采样。 一个示例是从具有以下密度的分布中采样随机量 $X$
$$\frac{1}{C}\cdot \frac{1}{1+|x-2{|}^3}$$
其中 $C={\int }_{-\infty }^{+\infty }\frac{1}{1+|x-2{|}^3}dx$，无法显式计算出来。

编写伪代码，通过接受-拒绝方法从分布（1.3.1）中采样随机量（提示：取 $M=5$ ）。 编写一个程序来调整接受-拒绝过程 1000 次，计算创建了多少个随机量，并绘制直方图。

import numpy as np
import matplotlib.pyplot as plt# Define the target distribution function f(x)
def target_distribution(x):return 1 / (1 + np.abs(x - 2) ** 3)# Define the auxiliary distribution (Cauchy distribution)
def auxiliary_distribution(x, x0, gamma):return 1 / (np.pi * gamma * (1 + ((x - x0) / gamma) ** 2))# Set the number of samples to generate
num_samples = 1000# Initialize arrays to store accepted and rejected samples
accepted_samples = []
rejected_samples = []# Parameters for the Cauchy distribution
x0 = 2  # Center of the Cauchy distribution
gamma = 1  # Scale parameter of the Cauchy distribution# Set the maximum value for the acceptance ratio
max_f_x = 5# Acceptance rate
acceptance_rate = 0# Generate samples using accept-reject method
for _ in range(num_samples):# Generate a random sample from the Cauchy distributionx_sample = np.random.cauchy(x0, gamma)# Generate a uniform random number for acceptance/rejectionu = np.random.uniform(0, max_f_x)# Calculate the acceptance probabilityacceptance_prob = target_distribution(x_sample) / (auxiliary_distribution(x_sample, x0, gamma))# Accept or reject the sampleif u < acceptance_prob:accepted_samples.append(x_sample)acceptance_rate += 1else:rejected_samples.append(x_sample)# Calculate the acceptance rate
acceptance_rate /= num_samples# Plot the histogram of accepted samples
plt.hist(accepted_samples, bins=30, density=True, alpha=0.5, label='Accepted Samples')
plt.plot(x_range, target_distribution(x_range), 'r', label='Target Distribution')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.title(f'Acceptance Rate: {acceptance_rate:.2%}')
plt.show()print(f'Number of accepted samples: {len(accepted_samples)}')