用MCMC采样的MATLAB代码？

以下是一个简单的 Metropolis-Hastings MCMC 采样的 MATLAB 代码示例：

% 目标分布的概率密度函数
function p = target_distribution(x)
    p = exp(-0.5*(x-1)^2) + exp(-0.5*(x+1)^2);
end

% MCMC采样
num_samples = 10000; % 采样数量
proposal_std = 1; % 建议分布的标准差
samples = zeros(1, num_samples); % 存储采样结果的数组
x = 0; % 初始化采样状态
for i = 1:num_samples
    % 从建议分布中生成建议状态
    x_proposed = x + proposal_std*randn();
    % 计算接受率
    acceptance_prob = min(1, target_distribution(x_proposed)/target_distribution(x));
    % 决定是否接受建议状态
    if rand() < acceptance_prob
        x = x_proposed;
    end
    samples(i) = x; % 将当前状态存储到采样结果数组中
end

% 绘制采样结果的直方图
histogram(samples, 'Normalization', 'pdf');
hold on;
% 绘制目标分布的概率密度函数
x_vals = linspace(-4, 4, 100);
y_vals = arrayfun(@target_distribution, x_vals);
plot(x_vals, y_vals, 'LineWidth', 2);
xlabel('x');
ylabel('Probability density');
legend('Sampled distribution', 'Target distribution');

这个代码的目标分布是一个标准正态分布的混合物，即 $p(x) = \frac{1}{2}\mathcal{N}(x | -1, 1) + \frac{1}{2}\mathcal{N}(x | 1, 1)$。在每次迭代中，它从一个高斯分布（即建议分布）中生成一个建议状态，然后计算接受率并决定是否接受建议状态。最终，它将采样结果存储在 `samples` 数组中，并绘制采样结果的直方图。