matlab代码实现MCMC波达角估计

以下是一个简单的 MATLAB 代码实现 MCMC 波达角估计：

% 定义观测数据
n = 100;
theta_true = 30;
sigma = 5;
x = randn(1, n) * sigma + theta_true;

% 定义先验分布
prior_mean = 0;
prior_sd = 20;

% 定义 MCMC 参数
niter = 10000;
burnin = 1000;

% 初始化 MCMC
theta = zeros(1, niter);
theta(1) = 0;
accept = 0;

% 进行 MCMC
for i = 2:niter
    % 从提议分布中抽样
    theta_prop = theta(i-1) + randn * 2;
    
    % 计算接受率
    prior_ratio = normpdf(theta_prop, prior_mean, prior_sd) / normpdf(theta(i-1), prior_mean, prior_sd);
    likelihood_ratio = exp(-0.5 * sum((x - theta_prop).^2) / sigma^2) / exp(-0.5 * sum((x - theta(i-1)).^2) / sigma^2);
    ratio = prior_ratio * likelihood_ratio;
    alpha = min(1, ratio);
    
    % 决定是否接受提议
    u = rand;
    if u <= alpha
        theta(i) = theta_prop;
        accept = accept + 1;
    else
        theta(i) = theta(i-1);
    end
end

% 去掉燃烧期
theta = theta(burnin+1:end);

% 输出结果
fprintf('Posterior mean: %.2f\n', mean(theta));
fprintf('Acceptance rate: %.2f%%\n', accept / (niter - 1) * 100);

这个代码实现了一个简单的 MCMC 链，其中提议分布是一个标准正态分布，先验分布是一个均值为 0，标准差为 20 的正态分布，似然函数是一个正态分布，方差为 5。运行代码后，将输出后验均值和接受率。