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

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

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

什么是频率学派？

在频率学派中，察看样本是随机的，而参数是固定的、未知的数量。

概率被解释为一个随机过程的许多观测的预期频率。

有一种想法是 "实在的"，例如，在预测鱼的生存环境时，盐度和温度之间的相互作用有一个回归系数？

什么是贝叶斯学派？

在贝叶斯办法中，概率被解释为对信念的主观掂量。

所有的变量--因变量、参数和假如都是随机变量。咱们用数据来确定一个预计的确定性（可信度）。

这种盐度X温度的相互作用反映的不是相对的，而是咱们对鱼的生存环境所理解的货色（实质上是粗率的）。

指标

频率学派

保障正确的误差概率，同时思考到抽样、样本大小和模型。

  • 毛病：须要对置信区间、第一类和第二类谬误进行简单的解释。
  • 长处：更具备外在的 "客观性 "和逻辑上的一致性。

贝叶斯学派

剖析更多的信息能在多大程度上进步咱们对一个零碎的意识。

  • 毛病：这都是对于信奉的问题! ...有重大影响。
  • 长处: 更直观的解释和施行，例如，这是这个假如的概率，这是这个参数等于这个值的概率。可能更靠近于人类天然地解释世界的形式。

理论利用中：为什么用贝叶斯

  • 具备无限数据的简单模型，例如层次模型，其中

  • 理论的先验常识非常少

贝叶斯法令：

一些典型的贝叶斯速记法。

留神:

  • 贝叶斯的最大问题在于确定先验散布。先验应该是什么？它有什么影响？

指标:

计算参数的后验散布：（|X）。

点估计是后验的平均值。

一个可信的区间是

你能够把它解释为一个参数在这个区间内的概率 。

计算

皮埃尔-西蒙-拉普拉斯（1749-1827）（见：Sharon Bertsch McGrayne: The Theory That Would Not Die)

  • 有些问题是可剖析的，例如二项式似然-贝塔先验。

    • 如果你有几个参数，而且是奇数散布，你能够用数值乘以/整合先验和似然（又称网格近似）。

      • 但如果你有很多参数，这是不可能实现的操作
  • 只管该实践能够追溯到1700年，甚至它对推理的解释也能够追溯到19世纪初，但它始终难以更宽泛地施行，直到马尔科夫链蒙特卡洛技术的倒退。

MCMC

MCMC的思维是对参数值i进行 "抽样"。

回顾一下，马尔科夫链是一个随机过程，它只取决于它的前一个状态，而且（如果是遍历的），会生成一个安稳的散布。

技巧 "是找到渐进地靠近正确散布的抽样规定（MCMC算法）。\
 

有几种这样的（相干）算法。

  • Metropolis-Hastings抽样
  • Gibbs 抽样
  • No U-Turn Sampling (NUTS)
  • Reversible Jump

一个一直倒退的文献和工作体系!

Metropolis-Hastings 算法

  1. 开始:
  2. 跳到一个新的候选地位:
  3. 计算后验:
  4. 如果
  5. 如果
  6. 转到第2步

Metropolis-Hastings: 硬币例子

你抛出了5个侧面。你对的最后 "猜想 "是

MCMC:

 p.old <- prior *likelihood while(length(thetas) <= n){  theta.new <- theta + rnorm(1,0,0.05)  p.new <- prior *likelihood   if(p.new > p.old | runif(1) < p.new/p.old){    theta <- theta.new    p.old <- p.new  }

画图:

hist(thetas[-(1:100)] )curve(6*x^5 )

采样链：调整、细化、多链

  • 那个 "朝向 "安稳的初始过渡被称为 "预烧期"，必须加以修整。

    • 怎么做？用眼睛看

  • 采样过程（显然）是自相干的。

    • 如何做？通常是用眼看，用acf()作为领导。\
       
  • 为了保障你收敛到正确的散布，你通常会从不同的地位取得多条链（例如4条）。
  • 无效样本量

MCMC 诊断法

R软件包帮忙剖析MCMC链。一个例子是线性回归的贝叶斯拟合（,,

plot(line)

预烧局部:

plot(line[[1]], start=10)

MCMC诊断法

查看后验散布（同时评估收敛性）。

density(line)

参数之间的关联性，以及链内的自相干关系

levelplot(line[[2]])acfplot(line)

统计摘要

运行MCMC的工具（在R外部）

逻辑Logistic回归：婴儿出世体重低

logitmcmc(low~age+as.factor(race)+smoke )

plot(mcmc)

MCMC与GLM逻辑回归的比拟

MCMC与GLM逻辑回归的比拟

对于这个利用，没有很好的理由应用贝叶斯建模，除非--你是 "贝叶斯主义者"。 你有对于回归系数的真正先验信息（这基本上是不太可能的）。

一个次要的毛病是 先验散布辣手的调整参数。

然而，MCMC能够拟合的一些更简单的模型（例如，档次的logit MCMChlogit）。

Metropolis-Hastings

Metropolis-Hastings很好，很简略，很广泛。然而对循环次数很敏感。而且可能太慢，因为它最终会回绝大量的循环。

Gibbs 采样

在Gibbs吉布斯抽样中，你不是用适当的概率承受/回绝，而是用适当的条件概率在参数空间中前进。 并从该散布中抽取一次。

而后你从新的条件散布中抽取下一个参数。

比Metropolis-Hastings快得多。无效样本量要高得多!

BUGS（OpenBUGS，WinBUGS）是应用吉布斯采样器的贝叶斯推理。

JAGS是 "吉布斯采样器"

其余采样器

汉密尔顿蒙特卡洛（HMC）--是一种梯度的Metropolis-Hastings，因而速度更快，对参数之间的关联性更好。

No-U Turn Sampler（NUTS）--因为不须要固定的长度，它的速度更快。这是STAN应用的办法（见http\://arxiv.org/pdf/1111.4246v1.pdf）。

(Hoffman and Gelman 2011)

其余工具

你可能想创立你本人的模型，应用贝叶斯MC进行拟合，而不是依赖现有的模型。为此，有几个工具能够抉择。

  • BUGS / WinBUGS / OpenBUGS (Bayesian inference Using Gibbs Sampling) - 贝叶斯抽样工具的鼻祖（自1989年起）。WinBUGS是专有的。OpenBUGS的支持率很低。
  • JAGS（Just Another Gibbs Sampler）承受一个用相似于R语言的语法编写的模型字符串，并应用吉布斯抽样从这个模型中编译和生成MCMC样本。能够在R中应用rjags包。
  • Stan（以Stanislaw Ulam命名）是一个相似于JAGS的相当新的程序--速度更快，更弱小，倒退迅速。从伪R/C语法生成C++代码。装置：http\://mc-stan.org/rstan.html**
  • Laplace’s Demon 所有的贝叶斯工具都在R中： http\://www.bayesian-inference.com/software

STAN

\
 

要用STAN拟合一个模型，步骤是:

  1. 为模型生成一个STAN语法伪代码（在JAGS和BUGS中雷同
  2. 运行一个R命令，用C++语言编译该模型
  3. 应用生成的函数来拟合你的数据

STAN示例--线性回归

STAN代码是R（例如，具备散布函数）和C（即你必须申明你的变量）之间的一种混合。每个模型定义都有三个块。

1.数据块:

  int n; //  vector[n] y; // Y 向量

这指定了你要输出的原始数据。在本例中，只有Y和X，它们都是长度为n的（数字）向量，是一个不能小于0的整数。

2. 参数块

  real beta1;  // slope

这些列出了你要预计的参数：截距、斜率和方差。

3. 模型块

    sigma ~ inv_gamma(0.001, 0.001);     yhat[i] <- beta0 + beta1 * (x[i] - mean(x));}    y ~ normal(yhat, sigma);

留神:

  • 你能够矢量化，但循环也同样快
  • 有许多散布（和 "平均值 "等函数）可用
请常常参阅手册！ https://github.com/stan-dev/stan/releases/download/v2.9.0/stan-reference-2.9.0.pdf

2. 在R中编译模型

你把你的模型保留在一个独自的文件中， 而后用stan_model()命令编译这个模型。

这个命令是把你形容的模型，用C++编码和编译一个NUTS采样器。置信我，本人编写C++代码是一件十分十分苦楚的事件（如果没有很多教训的话），而且它保障比R中的等同代码快得多。

留神：这一步可能会很慢。

3. 在R中运行该模型

这里的要害函数是sampling()。还要留神的是，为了给你的模型提供数据，它必须是列表的模式

模仿一些数据。

X <- runif(100,0,20)Y <- rnorm(100, beta0+beta1*X, sigma)

进行取样!

sampling(stan, Data)

这里有大量的输入，因为它计算了

print(fit, digits = 2)

MCMC诊断法

为了利用coda系列的诊断工具，你须要从STAN拟合对象中提取链，并将其从新创立为mcmc.list。

extract(stan.fitalply(chains, 2, mcmc)

最受欢迎的见解

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模型实现