乐趣区

关于数据挖掘:R语言贝叶斯MCMCGLM逻辑回归Rstan线性回归Metropolis-Hastings与Gibbs采样算法附代码数据

jiezi

原文链接：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.fit
alply(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 模型实现

退出移动版