关于数据挖掘:拓端tecdatR语言贝叶斯MetropolisHastings-Gibbs-吉布斯采样器估计变点指数分布附代码数据

54次阅读

共计 1881 个字符,预计需要花费 5 分钟才能阅读完成。

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

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

最近咱们被客户要求撰写对于吉布斯采样器的钻研报告,包含一些图形和统计输入。
指数分布是泊松过程中事件之间工夫的概率分布,因而它用于预测到下一个事件的等待时间,例如,您须要在公共汽车站期待的工夫,直到下一班车到了。

在本文中,咱们将应用指数分布,假如它的参数 λ,即事件之间的均匀工夫,在某个工夫点 k 产生了变动,即:

 

咱们的次要指标是应用 Gibbs 采样器在给定来自该散布的 n 个观测样本的状况下预计参数 λ、α 和 k。

吉布斯 Gibbs 采样器

Gibbs 采样器是 Metropolis-Hastings 采样器的一个特例,通常在指标是多元散布时应用。应用这种办法,链是通过从指标散布的边缘散布中采样生成的,因而每个候选点都被承受。

Gibbs 采样器生成马尔可夫链如下:

  1. 让  是 Rd 中的随机向量,在工夫 t=0 初始化 X(0)。
  2. 对于每次迭代 t=1,2,3,… 反复:
  • 设置 x1=X1(t-1)。
  • 对于每个 j=1,…,d:

    • 生成 X∗j(t) 从,其中  是给定 X(-j) 的 Xj 的单变量条件密度。
    • 更新 . 
  • 当每个候选点都被承受时,设置 . 
  • 减少 t。

贝叶斯公式

变点问题的一个简略公式假如 f 和 g 已知密度:

 

 其中 k 未知且 k=1,2,…,n。让 Yi 为公交车达到公交车站之间通过的工夫(以分钟为单位)。假如变动点产生在第 k 分钟,即:

 

 当 Y=(Y1,Y2,…,Yn) 时,似然 L(Y|k)由下式给出:

 

假如具备独立先验的贝叶斯模型由下式给出:

数据和参数的联结散布为:

其中,

正如我之前提到的,Gibbs 采样器的实现须要从指标散布的边缘散布中采样,因而咱们须要找到 λ、α 和 k 的残缺条件散布。你怎么能这样做?简略来说,您必须从下面介绍的连贯散布中抉择仅依赖于感兴趣参数的项并疏忽其余项。

λ 的残缺条件散布由下式给出:

 

α 的残缺条件散布由下式给出:

k 的残缺条件散布由下式给出:

计算方法

在这里,您将学习如何应用应用 R 的 Gibbs 采样器来预计参数 λ、α 和 k。

数据

首先,咱们从具备变动点的下一个指数分布生成数据:

set.seed(98712)
y <- c(rexp(25, rate = 2), rexp(35, rate = 10))

思考到公交车站的状况,一开始公交车均匀每 2 分钟一班,但从工夫 i =26 开始,公交车开始均匀每 10 分钟一班到公交车站。

Gibbs 采样器的实现

首先,咱们须要初始化 k、λ 和 α。

n <- length(y) # 样本的察看值的数量
lci <- 10000 # 链的大小
aba <- alpha <- k <- numeric(lcan)
k[1] <- sample(1:n,

当初,对于算法的每次迭代,咱们须要生成 λ(t)、α(t) 和 k(t),如下所示(记住如果 k+1>n 没有变动点):

for (i in 2:lcan){kt <- k[i-1]
        # 生成 lambda
        lambda[i] <- rgamma
        # 生成 α
              # 产生 k   
        for (j in 1:n) {L[j] <- ((lambda[i] / alpha[i





# 删除链条上的前 9000 个值
bunIn <- 9000 

后果

在本节中,咱们将介绍 Gibbs 采样器生成的链及其参数 λ、α 和 k 的散布。参数的实在值用红线示意。

下表显示了参数的理论值和应用 Gibbs 采样器取得的估计值的平均值:

res <- c(mean(k[-(1:bun)]), mean(lmba[-(1:burn)]), mean(apa[-(1:buI)]))
resfil

论断

从后果中,咱们能够得出结论,应用 R 中的 Gibbs 采样器取得的具备变点的指数分布对参数 k、λ 和 α 的估计值的平均值靠近于参数的理论值,然而咱们冀望更好预计。这可能是因为抉择了链的初始值或抉择了 λ 和 α 的先验散布。


最受欢迎的见解

1. 应用 R 语言进行 METROPLIS-IN-GIBBS 采样和 MCMC 运行

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

3.R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样

4.R 语言 BUGS JAGS 贝叶斯剖析 马尔科夫链蒙特卡洛办法(MCMC)采样

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

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

7.R 语言用 Rcpp 减速 Metropolis-Hastings 抽样预计贝叶斯逻辑回归模型的参数

8.R 语言应用 Metropolis- Hasting 抽样算法进行逻辑回归

9. R 语言中基于混合数据抽样 (MIDAS) 回归的 HAR-RV 模型预测 GDP 增长回归的 HAR-RV 模型预测 GDP 增长 ”)

正文完
 0