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

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

指数分布是泊松过程中事件之间时间的概率分布，因此它用于预测到下一个事件的等待时间，例如，您需要在公共汽车站等待的时间，直到下一班车到了（点击文末“阅读原文”获取完整代码数据）。

相关视频

在本文中，我们将使用指数分布，假设它的参数 λ ，即事件之间的平均时间，在某个时间点 k 发生了变化，即：

我们的主要目标是使用 Gibbs 采样器在给定来自该分布的 n 个观测样本的情况下估计参数 λ、α 和 k。

吉布斯Gibbs 采样器

Gibbs 采样器是 Metropolis-Hastings 采样器的一个特例，通常在目标是多元分布时使用。使用这种方法，链是通过从目标分布的边缘分布中采样生成的，因此每个候选点都被接受。

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

• 让  是 Rd 中的随机向量，在时间 t=0 初始化 X(0)。

• 对于每次迭代 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分钟一班到公交车站。

---

点击标题查阅往期内容

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

左右滑动查看更多

01

02

03

04

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、λ 和 α 的估计值的平均值接近于参数的实际值，但是我们期望更好估计。这可能是由于选择了链的初始值或选择了 λ 和 α的先验分布。

---

点击文末“阅读原文”

获取全文完整资料。

本文选自《R语言贝叶斯METROPOLIS-HASTINGS GIBBS 吉布斯采样器估计变点指数分布分析泊松过程车站等待时间》。

点击标题查阅往期内容

R语言马尔可夫MCMC中的METROPOLIS HASTINGS，MH算法抽样（采样）法可视化实例

python贝叶斯随机过程：马尔可夫链Markov-Chain，MC和Metropolis-Hastings，MH采样算法可视化

Python贝叶斯推断Metropolis-Hastings（M-H）MCMC采样算法的实现

Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

Matlab用BUGS马尔可夫区制转换Markov switching随机波动率模型、序列蒙特卡罗SMC、M H采样分析时间序列

R语言RSTAN MCMC：NUTS采样算法用LASSO 构建贝叶斯线性回归模型分析职业声望数据

R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机波动率SV模型、粒子滤波、Metropolis Hasting采样时间序列分析

R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

R语言贝叶斯MCMC：用rstan建立线性回归模型分析汽车数据和可视化诊断

R语言贝叶斯MCMC：GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例

R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数

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

R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病

R语言中贝叶斯网络（BN）、动态贝叶斯网络、线性模型分析错颌畸形数据

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

Python贝叶斯回归分析住房负担能力数据集

R语言实现贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析

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

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

R语言和STAN,JAGS：用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言贝叶斯线性回归和多元线性回归构建工资预测模型

R语言贝叶斯推断与MCMC：实现Metropolis-Hastings 采样算法示例

R语言stan进行基于贝叶斯推断的回归模型

R语言中RStan贝叶斯层次模型分析示例

R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

WinBUGS对多元随机波动率模型：贝叶斯估计与模型比较

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

R语言贝叶斯推断与MCMC：实现Metropolis-Hastings 采样算法示例

R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

视频：R语言中的Stan概率编程MCMC采样的贝叶斯模型

R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计