原文链接：http://tecdat.cn/?p=3772

原文出处：拓端数据部落公众号

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

创立测试数据

第一步，咱们创立一些测试数据，用来拟合咱们的模型。咱们假如预测变量和因变量之间存在线性关系，所以咱们用线性模型并增加一些乐音。

trueA <- 5trueB <- 0trueSd <- 10sampleSize <- 31 # 创立独立的x值x <- (-(sampleSize-1)/2):((sampleSize-1)/2)# 依据ax + b + N（0，sd）创立因变量y <-  trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd) plot(x,y, main="Test Data")

图

定义统计模型

下一步是指定统计模型。咱们曾经晓得数据是用x和y之间的线性关系y = a * x + b和带有标准差sd的正态误差模型N（0，sd）创立的，所以让咱们应用雷同的模型进行拟合，看看如果咱们能够检索咱们的原始参数值。

从模型中导出似然函数

为了预计贝叶斯剖析中的参数，咱们须要导出咱们想要拟合的模型的似然函数。似然函数是咱们冀望察看到的数据以咱们所看到的模型的参数为条件产生的概率（密度）。因而，鉴于咱们的线性模型y = b + a*x + N(0,sd)将参数（a, b, sd）作为输出，咱们必须返回在这个模型下取得上述测试数据的概率（这听起来比较复杂，正如你在代码中看到的，咱们只是计算预测值y = b + a*x与察看到的y之间的差别，而后咱们必须查找这种偏差产生的概率密度（应用dnorm）。

likelihood <- function(param){    a = param[1]    b = param[2]    sd = param[3]         pred = a*x + b     sumll = sum(singlelikelihoods)     (sumll)  }  slopevalues <- function(x){return(likelihood(c(x, trueB, trueSd)))}

斜率参数对数似然曲线

作为阐明，代码的最初几行绘制了斜率参数a的一系列参数值的似然函数。

为什么咱们应用对数

您留神到后果是似然函数中概率的对数，这也是我对所有数据点的概率求和的起因（乘积的对数等于对数之和）。咱们为什么要做这个？ 因为很多小概率乘以的可能性很快就会变得十分小（比方10 ^ -34）。在某些阶段，计算机程序存在数字四舍五入的问题。 

定义先验

第二步，与贝叶斯统计中一样，咱们必须为每个参数指定先验散布。为了不便起见，我对所有三个参数应用了均匀分布和正态分布。 无信息先验通常是1 / sigma（如果你想理解起因，请看这里）。

#先验散布prior <- function(param){    a = param[1]    b = param[2]    sd = param[3]    aprior =  (a, min=0, max=10, log = T)    bprior = dnorm(b, sd = 5, log = T) }

**后验

先验和似然性的乘积是MCMC将要解决的理论数量。这个函数被称为后验 。同样，在这里咱们应用加总，因为咱们应用对数。

posterior <- function(param){   return ( (param) + prior(param))}

 

**MCMC

接下来是Metropolis-Hastings算法。该算法最常见的利用之一（如本例所示）是从贝叶斯统计中的后验密度中提取样本。然而，原则上，该算法可用于从任何可积函数中进行采样。因而，该算法的目标是在参数空间中跳转，然而以某种形式使得在某一点上的概率与咱们采样的函数成比例（这通常称为指标函数）。在咱们的例子中，这是下面定义的后验。

这是通过

  1. 从随机参数值开始
  2. 依据称为提议函数的某个概率密度，抉择靠近旧值的新参数值
  3. 以概率p（新）/ p（旧）跳到这个新点，其中p是指标函数，p> 1示意跳跃

当咱们运行这个算法时，它拜访的参数的散布会收敛到指标散布p。那么，让咱们在R中失去 ：

########Metropolis算法# ################ proposalfunction <- function(param){    return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))} run_metropolis_MCMC <- function(startvalue, iterations){      for (i in 1:iterations){                   if (runif(1) < probab){            chain[i+1,] = proposal        }else{            chain[i+1,] = chain[i,]        }    }    return(chain)}  chain = run_metropolis_MCMC(startvalue, 10000) burnIn = 5000acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))

应用后验的对数可能在开始时有点凌乱，特地是当您查看计算承受概率的行时（probab = exp（后验散布（提议散布） - 后验（链[i，]）） ）。要了解咱们为什么这样做，请留神p1 / p2 = exp [log（p1）-log（p2）]。

算法的第一步可能受初始值的偏差，因而通常被抛弃用于进一步剖析 。要看的一个乏味的输入是承受率： 承受规范回绝提议的频率是多少？承受率能够受提议函数的影响：通常，提议越靠近，承受率越大。然而，十分高的承受率通常是有益的：这意味着算法“停留”在同一点 。能够证实，20％到30％的承受率对于典型利用来说是最佳的 。

 

###概要####################### par(mfrow = c(2,3))hist( [-(1:burnIn),1],nclass=30, , main="Posterior of a", xlab="True value = red line" )abline(v = mean(chain[-(1:burnIn),1]))  #进行比拟：summary(lm(y~x))

所失去的图应该看起来像上面的图。你看到咱们检索到了或多或少用于创立数据的原始参数，你还看到咱们在最高后验值四周失去了肯定的区域，这些后验值也有一些数据，这相当于贝叶斯的置信区间。

图： 上排显示斜率（a）的后验预计，截距（b）和误差的标准偏差（sd）。下一行显示马尔可夫链参数值。

 

还有问题吗？请在上面留言！

最受欢迎的见解

1.matlab应用贝叶斯优化的深度学习

2.matlab贝叶斯隐马尔可夫hmm模型实现

3.R语言Gibbs抽样的贝叶斯简略线性回归仿真

4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

6.Python用PyMC3实现贝叶斯线性回归模型

7.R语言应用贝叶斯 层次模型进行空间数据分析

8.R语言随机搜寻变量抉择SSVS预计贝叶斯向量自回归（BVAR）模型

9.matlab贝叶斯隐马尔可夫hmm模型实现