对目标函数target_dist = @(x) x(1).*x(2)-1500用Gibbs生成 MCMC链的代码

以下是使用Gibbs生成MCMC链的代码，可以求解目标函数target_dist = @(x) x(1).*x(2)-1500的最大值。

% 定义目标函数
target_dist = @(x) x(1).*x(2)-1500;

% 初始化参数
x1 = 1;
x2 = 1;

% 定义迭代次数和burn-in期
n_iter = 50000;
burn_in = 10000;

% 定义标准差
sigma = [1, 1];

% 初始化MCMC链
chain = zeros(n_iter, 2);

% 使用Gibbs抽样生成MCMC链
for i = 1:n_iter
    % 抽取x1
    x1 = normrnd(chain(i, 2), sigma(1));
    
    % 抽取x2
    x2 = normrnd(chain(i, 1), sigma(2));
    
    % 计算接受率
    a = min(1, target_dist([x1, x2])/target_dist([chain(i, 1), chain(i, 2)]));
    
    % 接受/拒绝
    u = rand();
    if u < a
        chain(i+1, :) = [x1, x2];
    else
        chain(i+1, :) = chain(i, :);
    end
end

% 去除burn-in期
chain = chain(burn_in+1:end, :);

% 计算最大值和最优参数
[max_val, ind] = max(target_dist(chain));
opt_param = chain(ind, :);

% 绘制MCMC链的轨迹
figure;
plot(chain(:, 1), chain(:, 2), '.');
xlabel('x1');
ylabel('x2');
title('MCMC Chain');

% 绘制目标函数的等高线图和最优点
figure;
x1_range = linspace(min(chain(:, 1)), max(chain(:, 1)), 100);
x2_range = linspace(min(chain(:, 2)), max(chain(:, 2)), 100);
[X, Y] = meshgrid(x1_range, x2_range);
Z = target_dist([X(:), Y(:)]);
Z = reshape(Z, size(X));
contour(X, Y, Z, 30);
hold on;
plot(opt_param(1), opt_param(2), 'r*', 'MarkerSize', 10);
xlabel('x1');
ylabel('x2');
title('Target Distribution');

代码中首先定义了目标函数target_dist，并初始化了参数x1和x2。接着设置了迭代次数n_iter和burn-in期，以及标准差sigma。然后初始化MCMC链chain，并使用Gibbs抽样生成MCMC链。在每次迭代中，抽取x1和x2，计算接受率a，接受/拒绝新状态，更新MCMC链。最后去除burn-in期，计算最大值和最优参数，并绘制MCMC链的轨迹和目标函数的等高线图。