dTpdWou2D {sdetorus}

Approximation of the transition probability density of the WN diffusion in 2D

Description

Computation of the transition probability density (tpd) for a WN diffusion (with diagonal diffusion matrix)

Usage

dTpdWou2D(x, x0, t, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)

Arguments

x

a matrix of dimension c(n, 2) containing angles. They all must be in [\pi,\pi) so that the truncated wrapping by maxK windings is able to capture periodicity.

x0

a matrix of dimension c(n, 2) containing the starting angles. They all must be in [\pi,\pi). If all x0 are the same, a matrix of dimension c(1, 2) can be passed for better performance.

t

a scalar containing the times separating x and x0.

alpha

vector of length 3 parametrizing the A matrix as in alphaToA.

mu

a vector of length 2 giving the mean.

sigma

vector of length 2 containing the square root of the diagonal of \Sigma, the diffusion matrix.

rho

correlation coefficient of \Sigma.

maxK

maximum absolute value of the windings considered in the computation of the WN.

expTrc

truncation for exponential: exp(x) with x <= -expTrc is set to zero. Defaults to 30.

Details

The function checks for positive definiteness. If violated, it resets A such that solve(A) %*% Sigma is positive definite.

See Section 3.3 in García-Portugués et al. (2019) for details. See dTpdWou for the general case (less efficient for 1D).

Value

A vector of size n containing the tpd evaluated at x.

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

Examples

set.seed(3455267)
alpha <- c(2, 1, -1)
sigma <- c(1.5, 2)
rho <- 0.9
Sigma <- diag(sigma^2)
Sigma[1, 2] <- Sigma[2, 1] <- rho * prod(sigma)
A <- alphaToA(alpha = alpha, sigma = sigma, rho = rho)
solve(A) %*% Sigma
mu <- c(pi, 0)
x <- t(euler2D(x0 = matrix(c(0, 0), nrow = 1), A = A, mu = mu,
               sigma = sigma, N = 500, delta = 0.1)[1, , ])

sum(sapply(1:49, function(i) log(dTpdWou(x = matrix(x[i + 1, ], ncol = 2),
                                         x0 = x[i, ], t = 1.5, A = A,
                                         Sigma = Sigma, mu = mu))))

sum(log(dTpdWou2D(x = matrix(x[2:50, ], ncol = 2),
                  x0 = matrix(x[1:49, ], ncol = 2), t = 1.5, alpha = alpha,
                  mu = mu, sigma = sigma, rho = rho)))

lgrid <- 56
grid <- seq(-pi, pi, l = lgrid + 1)[-(lgrid + 1)]
image(grid, grid, matrix(dTpdWou(x = as.matrix(expand.grid(grid, grid)),
                                 x0 = c(0, 0), t = 0.5, A = A,
                                 Sigma = Sigma, mu = mu),
                         nrow = 56, ncol = 56), zlim = c(0, 0.25),
      main = "dTpdWou")
image(grid, grid, matrix(dTpdWou2D(x = as.matrix(expand.grid(grid, grid)),
                                   x0 = matrix(0, nrow = 56^2, ncol = 2),
                                   t = 0.5, alpha = alpha, sigma = sigma,
                                   mu = mu),
                         nrow = 56, ncol = 56), zlim = c(0, 0.25),
      main = "dTpdWou2D")

x <- seq(-pi, pi, l = 100)
t <- 1
image(x, x, matrix(dTpdWou2D(x = as.matrix(expand.grid(x, x)),
                             x0 = matrix(rep(0, 100 * 2), nrow = 100 * 100,
                                         ncol = 2),
                             t = t, alpha = alpha, mu = mu, sigma = sigma,
                             maxK = 2, expTrc = 30),
                             nrow = 100, ncol = 100),
      zlim = c(0, 0.25))
points(stepAheadWn2D(x0 = rbind(c(0, 0)), delta = t / 500,
                     A = alphaToA(alpha = alpha, sigma = sigma), mu = mu,
                     sigma = sigma, N = 500, M = 1000, maxK = 2,
                     expTrc = 30))