全文链接：http://tecdat.cn/?p=19664 

最近咱们被客户要求撰写对于MCMC的钻研报告，包含一些图形和统计输入。

MCMC是从简单概率模型中采样的通用技术。

  1. 蒙特卡洛
  2. 马尔可夫链
  3. 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算法蕴含以下内容：

  1. 初始化：随机抉择一个初始状态x；
  2. 依据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独立采样

  1. 从某种状态开始xt。代码中的x。
  2. 在代码中提出一个新的状态x'候选

  3. 计算“承受概率”

  4. 从[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<- 54000burnin<- 4000shape<- 2.5rate<-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，因而这些均为logdensity/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算法

该算法是从 后验密度中采样最常见的贝叶斯统计利用之一 。

下面定义的后验。

  1. 从随机参数值开始
  2. 依据某个候选函数的概率密度，抉择一个靠近旧值的新参数值
  3. 以概率p（new）/ p（old）跳到这个新点，其中p是指标函数，并且p> 1也意味着跳跃
  4. 请留神，咱们有一个 对称的跳跃/ 候选散布 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 = 5000accep = 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 ??1Residual standard error: 9.78 on 29 degrees of freedomMultiple 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采样用于回归的贝叶斯预计