遵义网站建设oadmin,工程机械网官网,北京网络营销北京,制作免费个人网站连续随机量的生成-接受拒绝重采样 1. 接受-拒绝重采样方法2. Python编程实现 1. 接受-拒绝重采样方法
重采样方法由两个步骤组成#xff0c;第一个步骤提供近似分布的随机量#xff0c;第二个步骤是校正机制。 在下文中#xff0c;我们用 f f f 表示目标分布#xff0c;用… 连续随机量的生成-接受拒绝重采样 1. 接受-拒绝重采样方法2. Python编程实现 1. 接受-拒绝重采样方法
重采样方法由两个步骤组成第一个步骤提供近似分布的随机量第二个步骤是校正机制。 在下文中我们用 f f f 表示目标分布用 g g g 表示辅助分布。
本课程我们只研究经典的接受-拒绝方法。
我们的目标是从概率密度为 f f f的分布中抽取一个随机量 X X X这可能很难抽取。 我们引入辅助概率分布 g g g并采取以下假设 (A) 假设存在一些 M 0 M0 M0 使得 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 XY XY 如果 U f ( Y ) M g ( Y ) U\frac{f(Y)}{M g(Y)} UMg(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 M5 M5 。 编写一个程序来调整接受-拒绝过程 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, bins30, densityTrue, alpha0.5, labelAccepted Samples)
plt.plot(x_range, target_distribution(x_range), r, labelTarget Distribution)
plt.xlabel(x)
plt.ylabel(Density)
plt.legend()
plt.title(fAcceptance Rate: {acceptance_rate:.2%})
plt.show()print(fNumber of accepted samples: {len(accepted_samples)})
