simulation - 在 R 中模拟 Ogata 的细化算法
问题描述
我正在尝试完全按照https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf中的算法 3 中给出的 Ogata 的细化算法,使用他们指定的参数来生成图 2。
这是逐字复制它的代码。
# Simulation of a Univariate Hawkes Poisson with
# Exponential Kernel γ(u) = αe^(−βu), on [0, T].
# Based on: https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf
library(tidyverse)
# Initialize parameters that remains the same for all realizations
mu <- 1.2
alpha <- 0.6
beta <- 0.8
T <- 10
# Initialize other variables that change through realizations
bigTau <- vector('numeric')
bigTau <- c(bigTau,0)
s <- 0
n <- 0
lambda_vec_accepted <- c(mu)
# Begin loop
while (s < T) {
# -------------------------------
# Compute lambda_bar
# -------------------------------
sum_over_big_Tau <- 0
for (i in c(1:length(bigTau))) {
sum_over_big_Tau <- sum_over_big_Tau + alpha*exp(-beta*(s-bigTau[i]))
}
lambda_bar <- mu + sum_over_big_Tau
u <- runif(1)
# so that w ∼ exponential(λ_bar)
w <- -log(u)/lambda_bar
# so that s is the next candidate point
s <- s+w
# Generate D ∼ uniform(0,1)
D <- runif(1)
# -------------------------------
# Compute lambda_s
# -------------------------------
sum_over_big_Tau <- 0
for (i in c(1:length(bigTau))) {
sum_over_big_Tau <- sum_over_big_Tau + alpha*exp(-beta*(s-bigTau[i]))
}
lambda_s <- mu + sum_over_big_Tau
# -------------------------------
# Accepting with prob. λ_s/λ_bar
# -------------------------------
if (D*lambda_bar <= lambda_s ) {
lambda_vec_accepted <- c(lambda_vec_accepted,lambda_s)
# updating the number of points accepted
n <- n+1
# naming it t_n
t_n <- s
# adding t_n to the ordered set bigTau
bigTau <- c(bigTau,t_n)
}
}
df<-data.frame(x=bigTau,y=lambda_vec_accepted)
ggplot(df) + geom_line(aes(x=bigTau,y=lambda_vec_accepted))
然而,我设法得到(运行多次)的 lambda 与时间的图是这样的,与他们在图 2 中得到的图(呈指数下降)相去甚远。
我不确定我在做什么错误。如果有人可以提供帮助,那就太好了。这是我的研究所需要的。我来自生物学,所以如果我做了一些愚蠢的事情,请原谅。谢谢。
解决方案
推荐阅读
- java - SonarQube:不应该直接存储或返回可变成员
- model-view-controller - 在 Windows Server 2012、IIS 8.5 上未找到视图或其主视图或没有视图引擎支持搜索位置错误
- php - csrf_token() 在 l5-swagger 中是空的,除了 GET 请求之外不能做任何请求
- python - Python:如何全屏显示并判断玩家人数
- java - Spring引导属性到对象
- c# - 类型不能作为非托管结构封送
- java - Google AdMob 没有在我的应用程序上显示任何广告印象?
- assembly - 8086 中通用寄存器的区别:[bx] 有效,[cx] 无效?
- javascript - 如何强制 iframe 视频在 Android 设备上旋转到横向模式?
- javascript - Extjs 构建无法解决依赖 Ext.app.Application