MCMC的MATLAB代码

MCMC（Markov Chain Monte Carlo）是一种常用的随机抽样方法，用于估计高维概率分布中的参数。在MATLAB中，有许多库和工具箱支持MCMC算法，如`Statistics and Machine Learning Toolbox`。以下是一个简单的Metropolis-Hastings算法的示例代码，这是MCMC中最基础的方法之一：

% 导入所需的库
if ~isToolboxAvailable('Statistics')
    error('Please install Statistics and Machine Learning Toolbox for MCMC examples.');
end

% 定义目标密度函数 (假设我们有一个正态分布)
function pdf = targetDensity(x, mu, sigma)
    pdf = 1 / (sigma * sqrt(2 * pi)) * exp(- mu)^2) / (2 * sigma^2));
end

% 初始化参数
mu_init = 0;   % 假设正态分布的均值
sigma_init = 1; % 假设标准差
numSamples = 10000; % 想要采样的样本数量
 burnInPeriod = floor(numSamples * 0.1); % 燃烧期丢弃

% 初始化链
chain = zeros(burnInPeriod + numSamples, 2);
x_current = mu_init;
chain(1, :) = [x_current; 0]; % 第一个点初始化为0

% MCMC循环
for i = 2:burnInPeriod + numSamples
    % 提取当前样本
    x_current = chain(i - 1, 1);
    
    % 生成提案
    proposal_mu = x_current + randn() * 0.5; % 均值附近随机移动
    proposal_sigma = sigma_init + randn() * 0.1; % 标准差随机变化
    
    % 计算目标密度和提案密度
    pdf_current = targetDensity(x_current, mu_init, sigma_init);
    pdf_proposal = targetDensity(proposal_mu, mu_init, proposal_sigma);
    
    % Metropolis-Hastings接受概率
    acceptance_ratio = min(1, pdf_proposal / pdf_current);
    
    % 随机决定是否接受提案
    if rand < acceptance_ratio
        chain(i, :) = [proposal_mu; proposal_sigma];
    else
        chain(i, :) = [x_current; proposal_sigma]; % 如果拒绝，保持原状态
    end
end

% 删除燃烧期并计算均值和标准差
posterior_samples = chain(burnInPeriod + 1:end, :);
[mean_post, std_post] = meanstd(posterior_samples(:, 1), 'omitnan');

% 打印结果
disp(['Mean: ', num2str(mean_post)])
disp(['Standard Deviation: ', num2str(std_post)])

% 相关问题--
1. MCMC在什么情况下最适合使用？
2. 如何调整Metropolis-Hastings算法的参数来优化采样效率？
3. 在实际应用中，如何确定MCMC的停止条件？

请注意，这只是一个基础示例，实际应用可能需要更复杂的适应性步骤、自适应尺度参数或其他改进方法。如果你有特定的概率分布或问题想要模拟，代码将有所不同。