R语言贝叶斯MetropolisHastings Gibbs 吉布斯采样器估计变点指数分布分析泊松过程车站等待时间

原文链接：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。

数据

首先，我们从具有变化点的下一个指数分布生成数据：

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

考虑到公交车站的情况，一开始公交车平均每2分钟一班，但从时间i=26开始，公交车开始平均每10分钟一班到公交车站。

Gibbs采样器的实现

首先，我们需要初始化 k、λ 和 α。

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

现在，对于算法的每次迭代，我们需要生成 λ(t)、α(t) 和 k(t)，如下所示（记住如果 k+1>n 没有变化点）：

1.  
   for (i in 2:lcan){
2.  
   kt <- k[i-1]
3.  
   # 生成lambda
4.  
   lambda[i] <- rgamma
5.  
   # 生成α
6.  
   # 产生k
7.  
   for (j in 1:n) {
8.  
   L[j] <- ((lambda[i] / alpha[i
9.  
    
10.  
     
11.  
     
12.  
     
13.  
     
14.  
    # 删除链条上的前9000个值
15.  
    bunIn <- 9000

结果

在本节中，我们将介绍 Gibbs 采样器生成的链及其参数 λ、α 和 k 的分布。参数的真实值用红线表示。

下表显示了参数的实际值和使用 Gibbs 采样器获得的估计值的平均值：

1.  
   res <- c(mean(k[-(1:bun)]), mean(lmba[-(1:burn)]), mean(apa[-(1:buI)]))
2.  
   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增长