dPsTpd {sdetorus}

Wrapped Euler and Shoji–Ozaki pseudo-transition probability densities

Description

Wrapped pseudo-transition probability densities.

Usage

dPsTpd(x, x0, t, method = c("E", "SO", "SO2"), b, jac.b, sigma2, b1, b2,
  circular = TRUE, maxK = 2, vmApprox = FALSE, twokpi = NULL, ...)

Arguments

x

a matrix of dimension c(n, p). If a vector is provided, is assumed that p = 1.

x0

a matrix of dimension c(n, p). If all x0 are the same, a matrix of dimension c(1, p) can be passed for better performance. If a vector is provided, is assumed that p = 1.

t

time step between x and x0.

method

a string for choosing "E" (Euler), "SO" (Shoji–Ozaki) or "SO2" (Shoji–Ozaki with Ito's expansion in the drift) method.

b

drift function. Must return a matrix of the same size as x.

jac.b

jacobian of the drift function.

sigma2

diagonal of the diffusion matrix (if univariate, this is the square of the diffusion coefficient). Must return an object of the same size as x.

b1

first derivative of the drift function (univariate). Must return a vector of the same length as x.

b2

second derivative of the drift function (univariate). Must return a vector of the same length as x.

circular

flag to indicate circular data.

maxK

maximum absolute winding number used if circular = TRUE.

vmApprox

flag to indicate von Mises approximation to wrapped normal. See
momentMatchWnVm and scoreMatchWnBvm.

twokpi

optional matrix of winding numbers to avoid its recomputation. See details.

...

additional parameters passed to b, b1, b2, jac.b and sigma2.

Details

See Section 3.2 in García-Portugués et al. (2019) for details. "SO2" implements Shoji and Ozai (1998)'s expansion with for p = 1. "SO" is the same expansion, for arbitrary p, but considering null second derivatives.

twokpi is repRow(2 * pi * c(-maxK:maxK), n = n) if p = 1 and
as.matrix(do.call(what = expand.grid, args = rep(list(2 * pi * c(-maxK:maxK)), p))) otherwise.

Value

Output from mleOptimWrapper.

References

García-Portugués, E., Sørensen, M., Mardia, K. V. and Hamelryck, T. (2019) Langevin diffusions on the torus: estimation and applications. Statistics and Computing, 29(2):1–22. doi:10.1007/s11222-017-9790-2

Shoji, I. and Ozaki, T. (1998) A statistical method of estimation and simulation for systems of stochastic differential equations. Biometrika, 85(1):240–243. doi:10.1093/biomet/85.1.240

Examples

# 1D
grid <- seq(-pi, pi, l = 501)[-501]
alpha <- 1
sigma <- 1
t <- 0.5
x0 <- pi/2
# manipulate::manipulate({

  # Drifts
  b <- function(x) driftWn1D(x = x, alpha = alpha, mu = 0, sigma = sigma)
  b1 <- function(x, h = 1e-4) {
    l <- length(x)
    res <- driftWn1D(x = c(x + h, x - h), alpha = alpha, mu = 0,
                     sigma = sigma)
    drop(res[1:l] - res[(l + 1):(2 * l)])/(2 * h)
  }
  b2 <- function(x, h = 1e-4) {
    l <- length(x)
    res <- driftWn1D(x = c(x + h, x, x - h), alpha = alpha, mu = 0,
                     sigma = sigma)
    drop(res[1:l] - 2 * res[(l + 1):(2 * l)] +
          res[(2 * l + 1):(3 * l)]) / (h^2)
  }

  # Squared diffusion
  sigma2 <- function(x) rep(sigma^2, length(x))

  # Plot
  plot(grid, dTpdPde1D(Mx = length(grid), x0 = x0, t = t, alpha = alpha,
                       mu = 0, sigma = sigma), type = "l",
       ylab = "Density", xlab = "", ylim = c(0, 0.75), lwd = 2)
  lines(grid, dTpdWou1D(x = grid, x0 = rep(x0, length(grid)), t = t,
                       alpha = alpha, mu = 0, sigma = sigma), col = 2)
  lines(grid, dPsTpd(x = grid, x0 = x0, t = t, method = "E", b = b,
                     b1 = b1, b2 = b2, sigma2 = sigma2), col = 3)
  lines(grid, dPsTpd(x = grid, x0 = x0, t = t, method = "SO", b = b,
                     b1 = b1, b2 = b2, sigma2 = sigma2), col = 4)
  lines(grid, dPsTpd(x = grid, x0 = x0, t = t, method = "SO2", b = b,
                     b1 = b1, b2 = b2, sigma2 = sigma2),
        col = 5)
  lines(grid, dPsTpd(x = grid, x0 = x0, t = t, method = "E", b = b,
                     b1 = b1, b2 = b2, sigma2 = sigma2, vmApprox = TRUE),
        col = 6)
  lines(grid, dPsTpd(x = grid, x0 = x0, t = t, method = "SO", b = b,
                     b1 = b1, b2 = b2, sigma2 = sigma2, vmApprox = TRUE),
        col = 7)
  lines(grid, dPsTpd(x = grid, x0 = x0, t = t, method = "SO2", b = b,
                     b1 = b1, b2 = b2, sigma2 = sigma2, vmApprox = TRUE),
        col = 8)
  legend("topright", legend = c("PDE", "WOU", "E", "SO1", "SO2", "EvM",
                                "SO1vM", "SO2vM"), lwd = 2, col = 1:8)

# }, x0 = manipulate::slider(-pi, pi, step = 0.1, initial = -pi),
# alpha = manipulate::slider(0.1, 5, step = 0.1, initial = 1),
# sigma = manipulate::slider(0.1, 5, step = 0.1, initial = 1),
# t = manipulate::slider(0.1, 5, step = 0.1, initial = 1))

# 2D
grid <- seq(-pi, pi, l = 76)[-76]
alpha1 <- 2
alpha2 <- 1
alpha3 <- 0.5
sig1 <- 1
sig2 <- 2
t <- 0.5
x01 <- pi/2
x02 <- -pi/2
# manipulate::manipulate({

  alpha <- c(alpha1, alpha2, alpha3)
  sigma <- c(sig1, sig2)
  x0 <- c(x01, x02)

  # Drifts
  b <- function(x) driftWn2D(x = x, A = alphaToA(alpha = alpha,
                                                 sigma = sigma),
                             mu = rep(0, 2), sigma = sigma)
  jac.b <- function(x, h = 1e-4) {
    l <- nrow(x)
    res <- driftWn2D(x = rbind(cbind(x[, 1] + h, x[, 2]),
                               cbind(x[, 1] - h, x[, 2]),
                               cbind(x[, 1], x[, 2] + h),
                               cbind(x[, 1], x[, 2] - h)),
                     A = alphaToA(alpha = alpha, sigma = sigma),
                     mu = rep(0, 2), sigma = sigma)
    cbind(res[1:l, ] - res[(l + 1):(2 * l), ],
          res[2 * l + 1:l, ] - res[2 * l + (l + 1):(2 * l), ]) / (2 * h)
  }

  # Squared diffusion
  sigma2 <- function(x) matrix(sigma^2, nrow = length(x) / 2L, ncol = 2)

  # Plot
  old_par <- par(mfrow = c(3, 2))
  plotSurface2D(grid, grid, z = dTpdPde2D(Mx = length(grid),
                                          My = length(grid), x0 = x0,
                                          t = t, alpha = alpha,
                                          mu = rep(0, 2), sigma = sigma),
                levels = seq(0, 1, l = 20), main = "Exact")
  plotSurface2D(grid, grid,
                f = function(x) drop(dTpdWou2D(x = x,
                                               x0 = repRow(x0, nrow(x)),
                                               t = t, alpha = alpha,
                                               mu = rep(0, 2),
                                               sigma = sigma)),
                levels = seq(0, 1, l = 20), fVect = TRUE, main = "WOU")
  plotSurface2D(grid, grid,
                f = function(x) dPsTpd(x = x, x0 = rbind(x0), t = t,
                                       method = "E", b = b, jac.b = jac.b,
                                       sigma2 = sigma2),
                levels = seq(0, 1, l = 20), fVect = TRUE, main = "E")
  plotSurface2D(grid, grid,
                f = function(x) dPsTpd(x = x, x0 = rbind(x0), t = t,
                                       method = "SO", b = b, jac.b = jac.b,
                                       sigma2 = sigma2),
                levels = seq(0, 1, l = 20), fVect = TRUE, main = "SO")
  plotSurface2D(grid, grid,
                f = function(x) dPsTpd(x = x, x0 = rbind(x0), t = t,
                                       method = "E", b = b, jac.b = jac.b,
                                       sigma2 = sigma2, vmApprox = TRUE),
                levels = seq(0, 1, l = 20), fVect = TRUE, main = "EvM")
  plotSurface2D(grid, grid,
                f = function(x) dPsTpd(x = x, x0 = rbind(x0), t = t,
                                       method = "SO", b = b, jac.b = jac.b,
                                       sigma2 = sigma2, vmApprox = TRUE),
                levels = seq(0, 1, l = 20), fVect = TRUE, main = "SOvM")
  par(old_par)

# }, x01 = manipulate::slider(-pi, pi, step = 0.1, initial = -pi),
# x02 = manipulate::slider(-pi, pi, step = 0.1, initial = -pi),
# alpha1 = manipulate::slider(0.1, 5, step = 0.1, initial = 1),
# alpha2 = manipulate::slider(0.1, 5, step = 0.1, initial = 1),
# alpha3 = manipulate::slider(-5, 5, step = 0.1, initial = 0),
# sig1 = manipulate::slider(0.1, 5, step = 0.1, initial = 1),
# sig2 = manipulate::slider(0.1, 5, step = 0.1, initial = 1),
# t = manipulate::slider(0.01, 5, step = 0.01, initial = 1))