共计 6276 个字符，预计需要花费 16 分钟才能阅读完成。
全文链接：http://tecdat.cn/?p=19664
最近咱们被客户要求撰写对于 Metropolis-Hastings 采样的钻研报告，包含一些图形和统计输入。
MCMC 是从简单概率模型中采样的通用技术。
- 蒙特卡洛
- 马尔可夫链
- Metropolis-Hastings 算法
问题
如果须要计算有简单后验 pdf p（θ| y）的随机变量 θ 的函数 f（θ）的平均值或期望值。
您可能须要计算后验概率分布 p（θ）的最大值。
解决期望值的一种办法是从 p（θ）绘制 N 个随机样本，当 N 足够大时，咱们能够通过以下公式迫近期望值或最大值
将雷同的策略利用于通过从 p（θ| y）采样并取样本集中的最大值来找到 argmaxp（θ| y）。
解决办法
1.1 间接模仿
1.2 逆 CDF
1.3 回绝 / 承受抽样
如果咱们不晓得准确 / 标准化的 pdf 或非常复杂，则 MCMC 会派上用场。
马尔可夫链
为了模仿马尔可夫链，咱们必须制订一个 过渡核 T（xi，xj）。过渡核是从状态 xi 迁徙到状态 xj 的概率。
马尔可夫链的收敛性意味着它具备安稳散布 π。马尔可夫链的 统计 散布 是安稳的, 那么它 意味着散布 不会随着工夫的推移而扭转。
Metropolis 算法
对于一个 Markov 链是安稳的。基本上示意
处于状态 x 并转换为状态 x ’ 的概率必须等于处于状态 x ’ 并转换为状态 x 的概率
或者
办法是将转换分为两个子步骤；候选和承受回绝。
令 q（x’| x）示意 候选密度，咱们能够应用概率 α（x’| x）来调整 q。
候选散布 Q（X’| X）是给定的候选 X 的状态 X ’ 的条件概率，
和 承受散布 α（x’| x）的条件概率承受候选的状态 X ’-X’。咱们设计了承受概率函数，以满足具体的均衡。
该 转移概率 能够写成：
插入上一个方程式，咱们有
Metropolis-Hastings 算法
A 的抉择遵循以下逻辑。
在 q 下从 x 到 x ’ 的转移太频繁了。因而，咱们应该抉择 α（x | x’）=1。然而，为了满足 粗疏安稳，咱们有
下一步是抉择满足上述条件的承受。Metropolis-Hastings 是一种常见的 抉择：
即，当接受度大于 1 时，咱们总是承受，而当接受度小于 1 时，咱们将相应地回绝。因而，Metropolis-Hastings 算法蕴含以下内容：
- 初始化：随机抉择一个初始状态 x；
- 依据 q（x’| x）随机抉择一个新状态 x ’;
3. 承受依据 α（x’| x）的状态。如果不承受，则不会进行转移，因而无需更新任何内容。否则，转移为 x ’；
4. 转移到 2，直到生成 T 状态；
5. 保留状态 x，执行 2。
原则上，咱们从散布 P（x）提取保留的状态，因为步骤 4 保障它们是不相干的。必须依据候选散布等不同因素来抉择 T 的值。重要的是，尚不分明应该应用哪种散布 q（x’| x）；必须针对以后的特定问题进行调整。
属性
Metropolis-Hastings 算法的一个乏味个性是它 仅取决于比率
是候选样本 x ’ 与先前样本 xt 之间的概率，
是两个方向（从 xt 到 x ’，反之亦然）的候选密度之比。如果候选密度对称，则等于 1。
马尔可夫链从任意初始值 x0 开始，并且算法运行屡次迭代，直到“初始状态”被“遗记”为止。这些被抛弃的样本称为预烧（burn-in）。其余的 x 可承受值集代表散布 P（x）中的样本
Metropolis 采样
一个简略的 Metropolis-Hastings 采样
让咱们看看从 伽玛散布 模仿任意形态和比例参数，应用具备 Metropolis-Hastings 采样算法。
上面给出了 Metropolis-Hastings 采样器的函数。该链初始化为零，并在每个阶段都倡议应用 N（a / b，a /（b * b））个候选对象。
基于正态分布且均值和方差雷同 gamma 的 Metropolis-Hastings 独立采样
- 从某种状态开始 xt。代码中的 x。
- 在代码中提出一个新的状态 x ’ 候选
-
计算“承受概率”
- 从[0,1] 得出一些均匀分布的随机数 u；如果 u <α 承受该点，则设置 xt + 1 = x’。否则，回绝它并设置 xt + 1 = xt。
MH 可视化
set.seed(123)
for (i in 2:n) {can <- rnorm(1, mu, sig)
aprob <- min(1, (dgamma(can, a, b)/dgamma(x,
a, b))/(dnorm(can, mu, sig)/dnorm(x,
mu, sig)))
u <- runif(1)
if (u < aprob)
x <- can
vec[i] <- x
画图
设置参数。
nrep<- 54000
burnin<- 4000
shape<- 2.5
rate<-2.6
批改图，仅蕴含预烧期后的链
vec=vec[-(1:burnin)]
#vec=vec[burnin:length(vec)]
<!—->
par(mfrow=c(2,1)) # 更改主框架，在一帧中有多少个图形
plot(ts(vec), xlab="Chain", ylab="Draws")
abline(h = mean(vec), lwd="2", col="red" )
点击题目查阅往期内容
Python 用 MCMC 马尔科夫链蒙特卡洛、回绝抽样和 Metropolis-Hastings 采样算法
左右滑动查看更多
01
02
03
04
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.007013 0.435600 0.724800 0.843300 1.133000 3.149000
<!—->
var(vec[-(1:burnin)])
<!—->
[1] 0.2976507
初始值
第一个样本 vec 是咱们链的初始 / 起始值。咱们能够更改它，以查看收敛是否产生了变动。
x <- 3*a/b
vec[1] <- x
抉择计划
如果候选密度与指标散布 P（x）的形态匹配，即 q（x’| xt）≈P（x’）q（x’|），则该算法成果最佳。xt）≈P（x’）。如果应用正态候选密度 q，则在预烧期间必须调整方差参数 σ2。
通常，这是通过计算承受率来实现的，承受率是在最初 N 个样本的窗口中承受的候选样本的比例。
如果 σ2 太大，则承受率将非常低，因为候选可能落在概率密度低得多的区域中，因而 a1 将十分小，且链将收敛得十分慢。
示例 2：回归的贝叶斯预计
Metropolis-Hastings 采样用于贝叶斯预计回归模型。
设定参数
DGP 和图
# 创立独立的 x 值，大概为零
x <- (-(Size-1)/2):((Size-1)/2)
# 依据 ax + b + N（0，sd）创立相干值
y <- trueA * x + trueB + rnorm(n=Size,mean=0,sd=trueSd)
正态分布拟然
pred = a*x + b
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
为什么应用对数
似然函数中概率的对数，这也是我求和所有数据点的概率（乘积的对数等于对数之和）的起因。
咱们为什么要做这个？强烈建议这样做，因为许多小概率相乘的概率会变得很小。在某个阶段，计算机程序会陷入数值四舍五入或下溢问题。
因而，当您编写概率时，请始终应用对数
示例：绘制斜率 a 的似然曲线
# 示例：绘制斜率 a 的似然曲线
plot (seq(3, 7, by=.05), slopelikelihoods , type="l")
先验散布
这三个参数的均匀分布和正态分布。
# 先验散布
# 更改优先级，log 为 True，因而这些均为 log
density/likelihood
aprior = dunif(a, min=0, max=10, log = T)
bprior = dnorm(b, sd = 2, log = T)
sdprior = dunif(sd, min=0, max=30, log = T)
后验
先验和概率的乘积是 MCMC 将要解决的理论量。此函数称为后验函数。同样，这里咱们应用和，因为咱们应用对数。
posterior <- function(param){return (likelihood(param) + prior(param))
}
Metropolis 算法
该算法是从 后验密度中采样最常见的贝叶斯统计利用之一。
下面定义的后验。
- 从随机参数值开始
- 依据某个候选函数的概率密度，抉择一个靠近旧值的新参数值
- 以概率 p（new）/ p（old）跳到这个新点，其中 p 是指标函数，并且 p > 1 也意味着跳跃
- 请留神，咱们有一个 对称的跳跃 / 候选散布 q（x’| x）。
标准差 σ 是固定的。
所以承受概率等于
######## Metropolis 算法 ################
for (i in 1:iterations){probab = exp(posterior(proposal) - posterior(chain[i,]))
if (runif(1) < probab){chain[i+1,] = proposal
}else{chain[i+1,] = chain[i,]
}
施行
（e）输入承受的值，并解释。
chain = metrMCMC(startvalue, 5500)
burnIn = 5000
accep = 1-mean(duplicated(chain[-(1:burnIn),]))
算法的第一步可能会因初始值而有偏差，因而通常会被抛弃来进行进一步剖析（预烧期）。令人感兴趣的输入是承受率：候选多久被算法承受回绝一次？候选函数会影响承受率：通常，候选越靠近，承受率就越大。然而，十分高的承受率通常是有益的：这意味着算法在同一点上“停留”，这导致对参数空间（混合）的解决不够现实。
咱们还能够更改初始值，以查看其是否更改后果 / 是否收敛。
startvalue = c(4,0,10)
小结
V1 V2 V3
Min. :4.068 Min. :-6.7072 Min. : 6.787
1st Qu.:4.913 1st Qu.:-2.6973 1st Qu.: 9.323
Median :5.052 Median :-1.7551 Median :10.178
Mean :5.052 Mean :-1.7377 Mean :10.385
3rd Qu.:5.193 3rd Qu.:-0.8134 3rd Qu.:11.166
Max. :5.989 Max. : 4.8425 Max. :19.223
<!—->
# 比拟:
summary(lm(y~x))
<!—->
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-22.259 -6.032 -1.718 6.955 19.892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.1756 1.7566 -1.808 0.081 .
x 5.0469 0.1964 25.697 <2e-16 ***
---
Signif. codes: 0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1
Residual standard error: 9.78 on 29 degrees of freedom
Multiple R-squared: 0.9579, Adjusted R-squared: 0.9565
F-statistic: 660.4 on 1 and 29 DF, p-value: < 2.2e-16
<!—->
summary(lm(y~x))$sigma
<!—->
[1] 9.780494
<!—->
coefficients(lm(y~x))[1]
<!—->
(Intercept)
-3.175555
<!—->
coefficients(lm(y~x))[2]
<!—->
x
5.046873
总结：
### 总结: #######################
par(mfrow = c(2,3))
hist(chain[-(1:burnIn),1],prob=TRUE,nclass=30,col="109"
abline(v = mean(chain[-(1:burnIn),1]), lwd="2")
点击文末 “浏览原文”
获取全文残缺代码数据资料。
本文选自《R 语言 MCMC:Metropolis-Hastings 采样用于回归的贝叶斯预计》。
点击题目查阅往期内容
Python 用 MCMC 马尔科夫链蒙特卡洛、回绝抽样和 Metropolis-Hastings 采样算法 \
R 语言贝叶斯 METROPOLIS-HASTINGS GIBBS 吉布斯采样器预计变点指数分布剖析泊松过程车站等待时间 \
R 语言马尔可夫 MCMC 中的 METROPOLIS HASTINGS，MH 算法抽样（采样）法可视化实例 \
python 贝叶斯随机过程：马尔可夫链 Markov-Chain，MC 和 Metropolis-Hastings，MH 采样算法可视化 \
Python 贝叶斯推断 Metropolis-Hastings（M-H）MCMC 采样算法的实现 \
Metropolis Hastings 采样和贝叶斯泊松回归 Poisson 模型 \
Matlab 用 BUGS 马尔可夫区制转换 Markov switching 随机稳定率模型、序列蒙特卡罗 SMC、M H 采样剖析工夫序列 R 语言 RSTAN MCMC：NUTS 采样算法用 LASSO 构建贝叶斯线性回归模型剖析职业声望数据 \
R 语言 BUGS 序列蒙特卡罗 SMC、马尔可夫转换随机稳定率 SV 模型、粒子滤波、Metropolis Hasting 采样工夫序列剖析 \
R 语言 Metropolis Hastings 采样和贝叶斯泊松回归 Poisson 模型 \
R 语言贝叶斯 MCMC：用 rstan 建设线性回归模型剖析汽车数据和可视化诊断 \
R 语言贝叶斯 MCMC：GLM 逻辑回归、Rstan 线性回归、Metropolis Hastings 与 Gibbs 采样算法实例 \
R 语言贝叶斯 Poisson 泊松 - 正态分布模型剖析职业足球比赛进球数 \
R 语言用 Rcpp 减速 Metropolis-Hastings 抽样预计贝叶斯逻辑回归模型的参数 \
R 语言逻辑回归、Naive Bayes 贝叶斯、决策树、随机森林算法预测心脏病 \
R 语言中贝叶斯网络（BN）、动静贝叶斯网络、线性模型剖析错颌畸形数据 \
R 语言中的 block Gibbs 吉布斯采样贝叶斯多元线性回归 \
Python 贝叶斯回归剖析住房累赘能力数据集 \
R 语言实现贝叶斯分位数回归、lasso 和自适应 lasso 贝叶斯分位数回归剖析 \
Python 用 PyMC3 实现贝叶斯线性回归模型 \
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型 \
R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析 \
R 语言和 STAN,JAGS：用 RSTAN,RJAG 建设贝叶斯多元线性回归预测选举数据 \
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研 \
R 语言贝叶斯线性回归和多元线性回归构建工资预测模型 \
R 语言贝叶斯推断与 MCMC：实现 Metropolis-Hastings 采样算法示例 \
R 语言 stan 进行基于贝叶斯推断的回归模型 \
R 语言中 RStan 贝叶斯层次模型剖析示例 \
R 语言应用 Metropolis-Hastings 采样算法自适应贝叶斯预计与可视化 \
R 语言随机搜寻变量抉择 SSVS 预计贝叶斯向量自回归（BVAR）模型 \
WinBUGS 对多元随机稳定率模型：贝叶斯预计与模型比拟 \
R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样 \
R 语言贝叶斯推断与 MCMC：实现 Metropolis-Hastings 采样算法示例 \
R 语言应用 Metropolis-Hastings 采样算法自适应贝叶斯预计与可视化 \
视频：R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型 \
R 语言 MCMC:Metropolis-Hastings 采样用于回归的贝叶斯预计
