帮我写一个基于马尔可夫链蒙特卡罗法的放射源定位的MTLAB程序

好的，下面是一个基于马尔可夫链蒙特卡罗（MCMC）方法的放射源定位的MATLAB程序：

%% 参数设置
N = 10000; % MCMC迭代次数
a = 0.5; % 放射源辐射率
sigma = 1; % 探测器测量误差标准差
c = 1; % 探测器响应常数

%% 生成模拟数据
x_true = [1, 2]; % 放射源位置
n_detectors = 3; % 探测器数量
detectors = [0, 0; 0, 4; 3, 0]; % 探测器位置
y = zeros(n_detectors, 1);
for i = 1:n_detectors
    y(i) = c * a / (4 * pi * norm(detectors(i, :) - x_true)) + sigma * randn(); % 探测器测量值
end

%% MCMC算法
x = [0, 0]; % 初始位置
log_posterior = zeros(N, 1);
for i = 1:N
    x_proposed = x + randn(1, 2); % 从当前位置生成一个新位置
    if x_proposed(1) < 0 || x_proposed(2) < 0 % 约束条件：位置不能小于零
        log_posterior(i) = -inf;
        continue;
    end
    % 计算新位置的后验概率
    log_likelihood_proposed = 0;
    for j = 1:n_detectors
        log_likelihood_proposed = log_likelihood_proposed - 0.5 * log(2 * pi * sigma^2) - 0.5 * (y(j) - c * a / (4 * pi * norm(detectors(j, :) - x_proposed)))^2 / sigma^2;
    end
    log_prior_proposed = -log(1 + norm(x_proposed)); % 先验概率：位置服从均匀分布
    log_posterior_proposed = log_likelihood_proposed + log_prior_proposed;
    % Metropolis-Hastings接受拒绝准则
    if log_posterior_proposed >= log_posterior(i) || log(rand()) < exp(log_posterior_proposed - log_posterior(i))
        x = x_proposed;
        log_posterior(i) = log_posterior_proposed;
    else
        log_posterior(i) = log_posterior(i - 1);
    end
end

%% 结果可视化
figure;
scatter(x_true(1), x_true(2), 100, 'r', 'filled'); % 真实位置
hold on;
scatter(detectors(:, 1), detectors(:, 2), 50, 'b', 'filled'); % 探测器位置
scatter(x(1), x(2), 100, 'g', 'filled'); % 估计位置
axis equal;
legend('True location', 'Detectors', 'Estimated location');
title('Radiation source localization');

该程序使用了Metropolis-Hastings算法实现MCMC采样，其中先验概率为位置服从均匀分布，后验概率为测量值的似然与位置的先验的乘积。程序中约束条件为位置不能小于零。运行程序后，将得到一个可视化图像，其中红色点表示真实位置，蓝色点表示探测器位置，绿色点表示估计位置。