# 5  简单线性模型

data {
int<lower=0> N;   // number of data items
int<lower=0> K;   // number of predictors
matrix[N, K] x;   // predictor matrix
vector[N] y;      // outcome vector
}
parameters {
real alpha;           // intercept
vector[K] beta;       // coefficients for predictors
real<lower=0> sigma;  // error scale
}
model {
// alpha ~ normal(0, 1);  // prior
// beta ~ normal(0, 1);   // prior
y ~ normal(x * beta + alpha, sigma);  // likelihood
}

library(cmdstanr)
# 编译模型
mod_state_x77 <- cmdstan_model(
stan_file = "code/state_x77.stan",
compile = TRUE,
cpp_options = list(stan_threads = TRUE)
)
# 准备数据
state_x77 <- data.frame(
x = state.center$x, y = state.center$y,
state_name = state.name,
state_abb = state.abb,
state_region = state.region,
state_division = state.division,
state.x77, check.names = FALSE
)
state_x77_d <- list(
N = nrow(state_x77), # 观测记录的条数
K = 2,               # 协变量个数
x = state_x77[, c("Income", "Murder")], # N x 2 矩阵
y = state_x77[, "Life Exp"]             # N 向量
)
nchains <- 4 # 4 条迭代链
# 给每条链设置不同的参数初始值
inits_data <- lapply(1:nchains, function(i) {
list(
alpha = runif(1, 0, 1),
beta = runif(2, 1, 10),
sigma = runif(1, 1, 10)
)
})
# 采样拟合模型
fit_state_x77 <- mod_state_x77$sample( data = state_x77_d, # 观测数据 init = inits_data, # 迭代初值 iter_warmup = 1000, # 每条链预处理迭代次数 iter_sampling = 2000, # 每条链总迭代次数 chains = nchains, # 马尔科夫链的数目 parallel_chains = 1, # 指定 CPU 核心数，可以给每条链分配一个 threads_per_chain = 1, # 每条链设置一个线程 show_messages = FALSE, # 不显示迭代的中间过程 refresh = 0, # 不显示采样的进度 seed = 20190425 # 设置随机数种子，不要使用 set.seed() 函数 ) 模型参数估计结果如下： 参数的 $$\alpha,\beta_1,\beta_2$$ 后验均值估计与普通线性模型的拟合结果非常一致。采样结果可以直接传递给 bayesplot，绘制参数迭代的轨迹图和后验分布图。 library(ggplot2) library(bayesplot) library(patchwork) # 参数的后验分布 p1 <- mcmc_hist(fit_state_x77$draws(c("alpha", "beta[1]", "beta[2]")),
facet_args = list(
labeller = ggplot2::label_parsed,
strip.position = "top",
ncol = 1
)
) + theme_classic()
# 参数的迭代轨迹
p2 <- mcmc_trace(fit_state_x77\$draws(c("alpha", "beta[1]", "beta[2]")),
facet_args = list(
labeller = ggplot2::label_parsed,
strip.position = "top",
ncol = 1
)
) + theme_classic() + theme(legend.title = element_blank())
# 绘图
p1 | p2