用贝叶斯MCMC抽样解线性方程组 matlab程序举例

以下是使用贝叶斯MCMC抽样解线性方程组的matlab程序示例：

%设置线性方程组
A = [1 2 3; 4 5 6; 7 8 9];
b = [3; 6; 9];

%初始化MCMC参数
n = size(A,1);
N = 100000; %采样次数
x0 = zeros(n,1); %MCMC起始点
sigma = 0.1; %MCMC步长
burnin = 10000; %MCMC燃烧期

%定义先验分布
mu = zeros(n,1); %均值
Sigma = eye(n); %协方差矩阵

%定义目标分布
logp = @(x) -0.5*(A*x-b)'*(A*x-b)-0.5*(x-mu)'*(Sigma\(x-mu)); %目标分布的对数密度函数

%运行MCMC
samples = mcmc_rw(x0,logp,N,sigma,burnin);

%计算后验分布的均值和标准差
mean_x = mean(samples)';
std_x = std(samples)';

%输出结果
fprintf('后验分布的均值：\n');
disp(mean_x);
fprintf('后验分布的标准差：\n');
disp(std_x);

其中，`mcmc_rw`是一个自定义的函数，用于实现基于随机游走的MCMC抽样。完整的程序代码如下：

function samples = mcmc_rw(x0,logp,N,sigma,burnin)
%使用随机游走MCMC抽样
%输入参数：
%x0：MCMC起始点
%logp：目标分布的对数密度函数
%N：采样次数
%sigma：MCMC步长
%burnin：燃烧期
%输出参数：
%samples：采样结果

n = length(x0);
samples = zeros(n,N);

%初始化
x = x0;
logp_x = logp(x);

%采样
for i=1:N+burnin
    %随机游走
    x_prop = x+sigma*randn(n,1);
    logp_x_prop = logp(x_prop);
    
    %Metropolis-Hastings接受/拒绝步骤
    alpha = exp(logp_x_prop-logp_x);
    if rand<alpha
        x = x_prop;
        logp_x = logp_x_prop;
    end
    
    %记录采样结果
    if i>burnin
        samples(:,i-burnin) = x;
    end
end

end