贝叶斯分析

先查看各命令展开的示例,再按需展开其他选项及 R 实现。改写代码前,请阅读适用条件及差异说明。

bayesmh

正态及逆Gamma先验的贝叶斯高斯回归

Stata

bayesmh y x, likelihood(normal({var})) prior({y:}, normal(0,100)) prior({var}, igamma(1,1)) rseed(135)

R 包:MCMCpack

MCMCpack::MCMCregress(y ~ x, data = data, burnin = 100, mcmc = 1000, b0 = 0, B0 = .01, c0 = 2, d0 = 2, seed = 135, verbose = 0)

与 Stata 有差异 固定流程;请明确调整示例列名。仅简单正态回归模型;采样器和初始化不同。短链仅演示语法,不能证明收敛。

bayespredict

高斯回归响应的后验均值

Stata

bayespredict pmean, mean

R 包:MCMCpack

drop(cbind(1, data$x) %*% colMeans(as.matrix(draws)[, c("(Intercept)", "x")]))

与 Stata 有差异 固定流程;请明确调整示例列名。预测条件均值,观察噪声积分后为零。

bayesreps

高斯后验预测分布的重复结局

Stata

bayesreps yrep*, nreps(10)

R 包:MCMCpack

{ set.seed(136); keep <- seq(1, nrow(draws), length.out = 10); B <- as.matrix(draws)[keep, ]; means <- cbind(1, data$x) %*% t(B[, c("(Intercept)", "x")]); means + matrix(rnorm(length(means), sd = rep(sqrt(B[, "sigma2"]), each = nrow(data))), nrow(data)) }

与 Stata 有差异 固定流程;请明确调整示例列名。包含后验参数不确定性和观察噪声;不写入预测文件。

bayesstats

MCMC有效样本量

Stata

bayesstats ess

R 包:coda

coda::effectiveSize(draws)

与 Stata 有差异 固定流程;请明确调整示例列名。谱密度ESS估计量不同于部分Stata自相关或分批计算。

后验均值标准差和分位数

Stata

bayesstats summary

R 包:coda

summary(draws)

与 Stata 有差异 固定流程;请明确调整示例列名。汇总给定后验抽样;蒙特卡洛误差取决于估计方法。

变换系数的后验汇总

Stata

bayesstats summary (effect:exp({y:x}))

R 包:coda

summary(exp(draws[, "x", drop = FALSE]))

与 Stata 有差异 固定流程;请明确调整示例列名。先变换抽样再汇总;exp(E[b])不等于E[exp(b)]。

90%最高后验密度区间

Stata

bayesstats summary {y:x}, hpd clevel(90)

R 包:coda

coda::HPDinterval(draws[, "x", drop = FALSE], prob = .9)

与 Stata 有差异 固定流程;请明确调整示例列名。HPD算法和多峰分布处理可能不同。

bayestest

斜率为正的后验概率

Stata

bayestest interval {y:x}, lower(0)

R 包:base

mean(draws[, "x"] > 0)

与 Stata 有差异 固定流程;请明确调整示例列名。区间假设的蒙特卡洛概率,并非频率学派p值。

效应位于实质区间的后验概率

Stata

bayestest interval {y:x}, lower(1.8) upper(2.2)

R 包:base

mean(draws[, "x"] > 1.8 & draws[, "x"] < 2.2)

与 Stata 有差异 固定流程;请明确调整示例列名。实质区间应在查看后验结果前确定。