0

我正在尝试完全按照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 中得到的图（呈指数下降）相去甚远。

有效的 XHTML

我不确定我在做什么错误。如果有人可以提供帮助，那就太好了。这是我的研究所需要的。我来自生物学，所以如果我做了一些愚蠢的事情，请原谅。谢谢。

4

0 回答 0