原文链接: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分钟一班到公交车站。
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增长")