拒绝采样小记

section 1

拒绝采样最关键的部分,搞个矩形、向矩形里投点等等,所做的一切都是为了获得一个密度曲线所围成区域的均匀分布。只要能获得这样一个在密度曲线下满足均匀分布的样本,我们就可以获得与该密度曲线相匹配的随机变量的采样样本。方法是,只需把每个蓝点的横坐标提取出来,这些横坐标所构成的样本就是我们的目标样本。
在这里插入图片描述在这里插入图片描述
step1:用一个矩形将这个密度曲线套起来,把密度曲线 框在一个矩形里。
step2:向矩形里随机投点10000次(虚值)。随机投点意味着在矩形这块区域内,这些点是满足均匀分布的。
step3:有的点落在了密度曲线下侧,有的点落在了密度曲线的上侧。只保留密度曲线下侧的点。
step4:把每个蓝点的横坐标提取出来,这些横坐标所构成的样本就是我们的目标样本。

section2

上述用了一个矩形,该矩形就是一个满足均匀分布的建议分布,建议分布只是获得目标密度函数曲线下均匀分布样本的一个辅助工具。采用均匀分布作为建议分布有时效率很低,为什么这么说?从上例就可以看出,均匀分布的好多点(那些绿点)都被剔除了,造成了一种浪费。可以选择一些其他曲线来把密度曲线框起来,效率会提高一点。
在这里插入图片描述这里目标密度函数曲线为 h ( x ) h(x) h(x),对应于下图中的蓝线;建议分布密度曲线为 g ( x ) g(x) g(x),我们把 g ( x ) g(x) g(x)乘上一个常数因子 c c c,然后 用 c g ( x ) cg(x) cg(x)这条曲线 将目标密度曲线框起来。

假定满足 g ( x ) g(x) g(x)的随机变量易于采样,拒绝采样的步骤如下:

  1. g ( x ) g(x) g(x)采到一个样本数据,记 x ⋆ x^{\star} x,把它作为一个建议
  2. 要不要接受这个建议 作为满足 h ( x ) h(x) h(x)分布的一个样本数据呢?我们定义一个接受概率: α = h ( x ⋆ ) c g ( x ⋆ ) \displaystyle\alpha = \frac{h(x^{\star})}{c g(x^{\star})} α=cg(x)h(x)

也就是说,我们 α \alpha α的概率 接受 x ⋆ x^{\star} x 作为 h ( x ) h(x) h(x)分布的一个样本数据

实际操作中,取一个 U ( 0 , 1 ) U(0, 1) U(0,1)的随机数 μ \mu μ,如果 μ < α \mu<\alpha μ<α,就接受 x ⋆ x^{\star} x作为 h ( x ) h(x) h(x)的一个样本数据;否则,把它舍弃掉,回到1.步继续循环。最终可以获得一个样本。

  • 文章开头是一下子抽取10000个点,到后来怎么成了一个个抽了呢?其实它们是对应的,把蓝点去掉的过程 就相当于做是否拒绝判断的过程。
  • 如果有密度曲线下的均匀分布样本,就可以得到与密度曲线相匹配的分布的一个样本。
  • 如果建议分布的形状和目标分布越接近,采样的效率就越高。

参考文献:https://gaolei786.github.io/statistics/reject.html

拒绝采样是蒙特卡洛方法中的一种采样技术,用于从不易直接采样的目标分布中抽取样本。在拒绝采样中,我们首先选择一个容易采样的建议分布,这个分布需要比目标分布更加“宽松”,即在目标分布的所有值上都有较高的概率值。然后,从建议分布中抽取样本,并根据一定的拒绝率来决定是否接受这个样本,以确保最终样本符合目标分布。 以下是一个简单的拒绝采样MATLAB代码示例,目标分布是指数分布,而建议分布是均匀分布: ```matlab % 参数设置 target_density = @(x) exp(-x); % 目标分布的概率密度函数 proposal_density = @(x) ones(size(x)); % 建议分布的概率密度函数,这里用均匀分布 proposal_support = [0, 10]; % 建议分布的支持范围 n = 1000; % 希望抽取的样本数量 % 拒绝采样 accepted_samples = []; while length(accepted_samples) < n % 从建议分布中抽取样本 proposal_samples = proposal_support(1) + (proposal_support(2) - proposal_support(1)) * rand(1, n); % 计算目标分布与建议分布的概率密度比值 ratio = target_density(proposal_samples) ./ proposal_density(proposal_samples); % 生成接受概率 accept_prob = min(ratio, 1); % 根据接受概率决定是否接受样本 accept = rand(1, n) < accept_prob; accepted_samples = [accepted_samples; proposal_samples(accept)]; end % 移除空行 accepted_samples = accepted_samples'; accepted_samples = accepted_samples(~isempty(accepted_samples), :); % 绘制结果 histogram(accepted_samples, 'Normalization', 'pdf'); hold on; f = @(x) target_density(x); x_values = linspace(min(accepted_samples), max(accepted_samples), 100); plot(x_values, f(x_values), 'r', 'LineWidth', 2); title('拒绝采样结果'); xlabel('样本值'); ylabel('概率密度'); legend('样本分布', '目标分布'); ``` 在这段代码中,我们首先定义了目标分布和建议分布的概率密度函数,然后通过循环从建议分布中抽取样本,并计算每个样本被接受的概率。如果样本被接受,就将其加入到最终样本集中。最后,我们绘制了样本的直方图和目标分布的曲线图,以可视化结果。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值