R/drawUniMMPPIntensity.R

In ppdiag: Diagnosis and Visualizations Tools for Temporal Point Processes

#' Draw the intensity of the Markov-modulated Poisson Process(MMPP)
#'
#' Take a mmpp object and draw its intensity accordingly
#'
#' @param mmpp a mmpp object including its transition probability matrix, lambda0, delta, and c.
#' @param simulation the simulated Markov-modulated Poisson Process(MMPP)
#' @param add logical; if TRUE add to an already existing plot;
#'  if NA start a new plot taking the defaults
#'  for the limits and log-scaling of the x-axis from the previous plot.
#'   Taken as FALSE (with a warning if a different value is supplied)
#'   if no graphics device is open.
#' @param color A specification for the default plotting color.
#' @param fit a boolean indicating whether to fit the events provided
#' @param int_title title of the plot.
#' @importFrom graphics plot
#' @importFrom graphics points
#' @importFrom graphics legend
#' @return no return value, intensity plot of Markov-modulated Poisson process
#' @export
#' @examples
#' Q <- matrix(c(-0.4, 0.4, 0.2, -0.2), ncol = 2, byrow = TRUE)
#' x <- pp_mmpp(Q, delta = c(1 / 3, 2 / 3), lambda0 = 0.9, c = 1.2)
#' y <- pp_simulate(x, n = 10)
#' drawUniMMPPIntensity(x, y)
drawUniMMPPIntensity <- function(mmpp, simulation,
                                 add = FALSE, color = 1, fit = FALSE,
                                 int_title = "Intensity Plot of MMPP") {
  # input mmpp: mmpp object generated by mmpp.R
  events <- simulation$events
  event_state <- simulation$zt
  state <- simulation$z
  state_time <- simulation$x
  lambda0 <- mmpp$lambda0
  c <- mmpp$c
  lambda1 <- lambda0 * (1 + c)

  n <- length(events)
  m <- length(state)

  yupper <- max(c(lambda0, lambda1)) + 1

  if (add == FALSE) {
    graphics::plot(0, 0,
      xlim = c(0, max(events)), ylim = c(0, yupper * 1.5), type = "n",
      xlab = "Time", ylab = "Intensity",
      main = int_title
    )
    ## should be related to state
    graphics::points(events[-1], rep(lambda0 / 2, n - 1),
      cex = 0.6,
      pch = ifelse(event_state[-1] == 1, 16, 1), col = "blue"
    )
    points(state_time, rep(lambda0, m), cex = 0.6, pch = 4, col = "red")
    for (s in state_time[2:(length(state_time) - 1)]) {
      graphics::segments(x0 = s, y0 = lambda0, y1 = lambda1, lty = 2)
    }
  }
  for (i in 1:(m - 1)) {
    if (state[i] == 2) {
      poisson_time <- events[events >= state_time[i] &
        events < state_time[i + 1]]
      if (i == 1) poisson_time <- poisson_time[-1]
      poisson_obj <- list(lambda = lambda0)
      class(poisson_obj) <- "hpp"
      drawHPPIntensity(poisson_obj,
        events = poisson_time, add = TRUE,
        start = state_time[i], end = state_time[i + 1], color = color, plot_events = FALSE
      )
    } else {
      poisson_time <- events[events >= state_time[i] &
        events < state_time[i + 1]]
      if (i == 1) poisson_time <- poisson_time[-1]
      poisson_obj <- list(lambda = lambda1)
      class(poisson_obj) <- "hpp"
      drawHPPIntensity(poisson_obj,
        events = poisson_time, add = TRUE,
        start = state_time[i], end = state_time[i + 1], color = color, plot_events = FALSE
      )
    }
  }
  if (add == FALSE) {
    graphics::legend("topleft", c(
      "Poisson state_1", "Poisson state_2",
      "state change point"
    ),
    col = c("blue", "blue", "red"),
    pch = c(16, 1, 4)
    )
  } else {
    legend("topright", c("True", "Estimation"),
      col = c("black", color),
      lty = c(1, 1)
    )
  }
}