连续随机量的生成-接受拒绝重采样
1. 接受-拒绝重采样方法
重采样方法由两个步骤组成,第一个步骤提供近似分布的随机量,第二个步骤是校正机制。 在下文中,我们用 f f f 表示目标分布,用 g g g 表示辅助分布。
本课程我们只研究经典的接受-拒绝方法。
我们的目标是从概率密度为
f
f
f的分布中抽取一个随机量
X
X
X,这可能很难抽取。 我们引入辅助概率分布
g
g
g并采取以下假设:
(A) 假设存在一些
M
>
0
M>0
M>0 使得
f
(
x
)
g
(
x
)
≤
M
对于所有
x
。
\frac{f(x)}{g(x)} \leq M \quad \text { 对于所有 } x \text {。 }
g(x)f(x)≤M 对于所有 x。
接受-拒绝方法的伪代码:
- 从分布 g g g 中采样随机数量 Y
- 从 U ( [ 0 , 1 ] ) U([0,1]) U([0,1]) 中抽取随机数量 U U U。
-
- 如果 U ⩽ f ( Y ) M g ( Y ) U \leqslant \frac{f(Y)}{M g(Y)} U⩽Mg(Y)f(Y),接受,即 X = Y X=Y X=Y
-
- 如果 U > f ( Y ) M g ( Y ) U>\frac{f(Y)}{M g(Y)} U>Mg(Y)f(Y),则拒绝并重复该过程。
该算法的理论证明:
让我们考虑一维情况。 对于任何
x
∈
R
x \in \mathbb{R}
x∈R,考虑
P
(
X
⩽
x
)
=
P
(
Y
⩽
x
∣
f
(
Y
)
M
g
(
Y
)
⩾
U
)
=
P
(
Y
⩽
x
,
f
(
Y
)
M
g
(
Y
)
⩾
U
)
P
(
f
(
Y
)
M
g
(
Y
)
⩽
U
)
=
∫
−
∞
x
(
∫
0
f
(
Y
)
M
g
(
y
)
d
u
)
g
(
y
)
d
y
∫
−
∞
+
∞
(
∫
0
f
(
y
)
M
g
(
y
)
d
u
)
g
(
y
)
d
y
(
∵
f
r
a
c
f
(
y
)
M
g
(
y
)
⩽
1
对于所有
y
)
=
∫
−
∞
x
f
(
y
)
d
y
∫
−
∞
+
∞
f
(
y
)
d
y
\begin{aligned} P(X\leqslant x) & =P\left(Y \leqslant x \mid \frac{f(Y)}{M g(Y)} \geqslant U\right) \\ & =\frac{P\left(Y \leqslant x, \frac{f(Y)}{M g(Y)} \geqslant U\right)}{P\left(\frac{f(Y)}{ M g(Y)} \leqslant U\right)} \\ & =\frac{\int_{-\infty}^x\left(\int_0^{\frac{f(Y)}{M g(y)}} d u\right) g(y) d y}{\int_ {-\infty}^{+\infty}\left(\int_0^{\frac{f(y)}{M g(y)}} d u\right) g(y) d y}\left(\because \ frac{f(y)}{M g(y)} \leqslant 1 \text { 对于所有 } y\right) \\ & =\frac{\int_{-\infty}^x f(y) d y}{\int_{-\infty}^{+\infty} f(y) d y} \end{aligned}
P(X⩽x)=P(Y⩽x∣Mg(Y)f(Y)⩾U)=P(Mg(Y)f(Y)⩽U)P(Y⩽x,Mg(Y)f(Y)⩾U)=∫−∞+∞(∫0Mg(y)f(y)du)g(y)dy∫−∞x(∫0Mg(y)f(Y)du)g(y)dy(∵ fracf(y)Mg(y)⩽1 对于所有 y)=∫−∞+∞f(y)dy∫−∞xf(y)dy
因此,
X
X
X 具有密度为
f
f
f 的分布。
2. Python编程实现
接受-预测方法可用于对分布乘上一个常数的随机量进行采样。 一个示例是从具有以下密度的分布中采样随机量
X
X
X
1
C
⋅
1
1
+
∣
x
−
2
∣
3
\frac{1}{C} \cdot \frac{1}{1+|x-2|^3}
C1⋅1+∣x−2∣31
其中
C
=
∫
−
∞
+
∞
1
1
+
∣
x
−
2
∣
3
d
x
C=\int_{-\infty}^{+\infty} \frac{1}{1+|x-2|^3} d x
C=∫−∞+∞1+∣x−2∣31dx,无法显式计算出来。
编写伪代码,通过接受-拒绝方法从分布(1.3.1)中采样随机量(提示:取 M = 5 M=5 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 auxiliary distribution (Cauchy distribution)
x_sample = np.random.standard_cauchy()
# Scale and shift the sample according to your Cauchy parameters
x_sample = x_sample * gamma + x0
# Generate a uniform random number for acceptance/rejection
u = np.random.uniform(0, max_f_x)
# Calculate the acceptance probability
acceptance_prob = target_distribution(x_sample) / (auxiliary_distribution(x_sample, x0, gamma))
# Accept or reject the sample
if u < acceptance_prob:
accepted_samples.append(x_sample)
acceptance_rate += 1
else:
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')
x_range = np.linspace(-5, 10, 1000)
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)}')