dTpdWou1D {sdetorus}

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

Description

Computation of the transition probability density (tpd) for a WN diffusion.

Usage

dTpdWou1D(x, x0, t, alpha, mu, sigma, maxK = 2L, expTrc = 30,
  vmApprox = 0L, kt = 0, logConstKt = 0)

Arguments

x

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

x0

a vector of length n containing the starting angles. They all must be in [\pi,\pi).

t

a scalar containing the times separating x and x0.

alpha

drift parameter.

mu

mean parameter. Must be in [\pi,\pi).

sigma

diffusion coefficient.

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.

vmApprox

whether to use the von Mises approximation to a wrapped normal (1) or not (0, default).

kt

concentration for the von Mises, a suitable output from momentMatchWnVm (see examples).

logConstKt

the logarithm of the von Mises normalizing constant associated to the concentration kt (see examples)

Details

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

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

t <- 0.5
alpha <- 1
mu <- 0
sigma <- 1
x0 <- 0.1
dTpdWou1D(x = seq(-pi, pi, l = 10), x0 = rep(x0, 10), t = t, alpha = alpha,
          mu = mu, sigma = sigma, vmApprox = 0)

# von Mises approximation
kt <- scoreMatchWnVm(sigma2 = sigma^2 * (1 - exp(-2 * alpha * t)) / (2 * alpha))
dTpdWou1D(x = seq(-pi, pi, l = 10), x0 = rep(x0, 10), t = t, alpha = alpha,
          mu = mu, sigma = sigma, vmApprox = 1, kt = kt,
          logConstKt = -log(2 * pi * besselI(x = kt, nu = 0,
                                             expon.scaled = TRUE)))