马尔科夫蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法是一种强大的工具,用于从复杂的概率分布中抽取样本。对于非常规的概率密度函数(PDF),MCMC方法尤其有用,因为这些分布可能难以直接采样。其中,Metropolis-Hastings算法是MCMC方法中最常用的一种。
Metropolis-Hastings算法的基本步骤
-
初始化:选择一个初始状态 \(x_0\)。
-
提议分布:选择一个提议分布 \(q(x'|x)\) ,用于生成候选样本 \(x'\) 。
-
接受率计算:计算接受率 \(\alpha\) :
\(\alpha = \min\left(1, \frac{p(x') q(x|x')}{p(x) q(x'|x)}\right)\)
其中 \(p(x)\) 是目标分布。
-
采样决策:生成一个均匀分布的随机数 \(u\) 。如果 \(u < \alpha\),则接受 \(x'\) 作为新的样本;否则,保持当前样本不变。
-
重复:重复上述步骤,直到获得足够多的样本。
MATLAB实现
使用Metropolis-Hastings算法从非常规概率密度函数中抽取样本的MATLAB示例。
定义目标分布
假设我们有一个非常规的概率密度函数 \(p(x)\),例如:
\(p(x) = \frac{1}{Z} \exp\left(-\frac{(x-3)^2}{2} - \frac{(x+3)^2}{2}\right)\)
其中 \(Z\) 是归一化常数。
function p = target_pdf(x)% 目标概率密度函数p = exp(-((x-3).^2)/2 - ((x+3).^2)/2);
end
Metropolis-Hastings算法实现
function [samples, acceptance_rate] = metropolis_hastings(target_pdf, num_samples, initial_state, proposal_std)% 输入参数:% target_pdf - 目标概率密度函数% num_samples - 需要生成的样本数量% initial_state - 初始状态% proposal_std - 提议分布的标准差% 初始化samples = zeros(1, num_samples);current_state = initial_state;accepted = 0;% 提议分布为正态分布proposal_dist = @(x) normpdf(x, 0, proposal_std);% Metropolis-Hastings算法for i = 1:num_samples% 生成候选样本candidate = current_state + proposal_std * randn;% 计算接受率alpha = min(1, (target_pdf(candidate) * proposal_dist(current_state - candidate)) / ...(target_pdf(current_state) * proposal_dist(candidate - current_state)));% 决定是否接受候选样本if rand < alphacurrent_state = candidate;accepted = accepted + 1;end% 保存样本samples(i) = current_state;end% 计算接受率acceptance_rate = accepted / num_samples;
end
示例
% 参数设置
num_samples = 10000; % 需要生成的样本数量
initial_state = 0; % 初始状态
proposal_std = 1; % 提议分布的标准差% 调用Metropolis-Hastings算法
[samples, acceptance_rate] = metropolis_hastings(@target_pdf, num_samples, initial_state, proposal_std);% 显示结果
figure;
histogram(samples, 100);
title('Sampled Distribution');
xlabel('x');
ylabel('Frequency');disp(['Acceptance Rate: ', num2str(acceptance_rate)]);
参考代码 使用马尔科夫蒙特卡洛方法对非常规的概率密度函数进行样本抽取 youwenfan.com/contentcnl/78939.html
说明
- 目标分布:目标分布 ( p(x) ) 可以是任意复杂的函数,只要能够计算其值即可。
- 提议分布:提议分布 ( q(x'|x) ) 通常选择为正态分布,但也可以根据问题选择其他分布。
- 接受率:接受率是一个重要的指标,通常希望接受率在0.2到0.5之间。如果接受率过高或过低,可以通过调整提议分布的标准差来优化。
上述代码,可以使用Metropolis-Hastings算法从非常规的概率密度函数中抽取样本,并通过直方图观察样本分布。