贝叶斯mcmc matlab实现

以下是一个简单的贝叶斯MCMC的Matlab实现：

假设我们要对以下模型进行贝叶斯推断：

y ~ N(θ,σ²)

θ ~ N(0,τ²)

σ² ~ IG(a,b)

其中，N表示正态分布，IG表示逆Gamma分布，a和b是逆Gamma分布的参数。

我们的目标是估计参数θ和σ²的后验分布。

首先，我们需要编写一个函数来计算我们的后验概率密度函数（PDF）：

function [log_posterior] = log_posterior(theta, sigma, y, tau, a, b)
log_prior = log(normpdf(theta, 0, tau)) + log(invgampdf(sigma, a, b));
log_likelihood = log(normpdf(y, theta, sigma^2));
log_posterior = sum(log_prior) + sum(log_likelihood);
end

其中，log_prior计算先验概率密度函数，log_likelihood计算似然函数，log_posterior计算后验概率密度函数。

接下来，我们将使用Metropolis-Hastings算法来生成参数的后验分布样本。我们将从一个初始值开始，然后使用下面的步骤来生成样本：

1. 从一个建议分布中抽取一个新的参数值（在这里，我们将使用正态分布作为建议分布）。

2. 计算接受概率（接受率）：

a. 计算新参数值的后验概率密度函数和旧参数值的后验概率密度函数。

b. 将两个概率密度函数相除，取对数。

c. 取此对数值和1的最小值，作为接受率。

3. 从均匀分布中抽取一个随机数，如果小于接受率，则接受新参数值。否则，保留旧参数值。

4. 重复步骤1到3，直到生成足够数量的样本。

下面是Matlab代码：

% 生成数据
n = 100;
theta_true = 2;
sigma_true = 1;
y = normrnd(theta_true, sigma_true, n, 1);

% 设置MCMC参数
num_samples = 10000;
burn_in = 1000;
thin = 10;
proposal_sd = [0.1, 0.1];

% 设置先验分布参数
tau = 1;
a = 1;
b = 1;

% 初始化参数
theta = zeros(num_samples, 1);
sigma = zeros(num_samples, 1);
theta(1) = 0;
sigma(1) = 1;

% 运行MCMC
for i = 2:num_samples
% 提议新的参数值
theta_prop = normrnd(theta(i-1), proposal_sd(1));
sigma_prop = normrnd(sigma(i-1), proposal_sd(2));

% 计算接受率
log_posterior_prop = log_posterior(theta_prop, sigma_prop, y, tau, a, b);
log_posterior_old = log_posterior(theta(i-1), sigma(i-1), y, tau, a, b);
log_acceptance_prob = log_posterior_prop - log_posterior_old;
acceptance_prob = min(1, exp(log_acceptance_prob));

% 决定是否接受新的参数值
if rand() < acceptance_prob
theta(i) = theta_prop;
sigma(i) = sigma_prop;
else
theta(i) = theta(i-1);
sigma(i) = sigma(i-1);
end
end

% 去除burn-in期和thin样本
theta = theta(burn_in+1:thin:end);
sigma = sigma(burn_in+1:thin:end);

% 打印参数估计
disp(['Mean of theta: ', num2str(mean(theta))]);
disp(['Mean of sigma: ', num2str(mean(sigma))]);
disp(['95% credible interval of theta: [', num2str(prctile(theta, 2.5)), ',', num2str(prctile(theta, 97.5)), ']']);
disp(['95% credible interval of sigma: [', num2str(prctile(sigma, 2.5)), ',', num2str(prctile(sigma, 97.5)), ']']);

注意，我们对生成的样本进行了去除burn-in期和thin样本的处理，以便得到更准确的参数估计。

这是一个简单的贝叶斯MCMC的Matlab实现，可能需要更复杂的模型和更复杂的算法来处理更具挑战性的问题。