| Title: | Statistical Tools for Toroidal Diffusions |
|---|---|
| Description: | Implementation of statistical methods for the estimation of toroidal diffusions. Several diffusive models are provided, most of them belonging to the Langevin family of diffusions on the torus. Specifically, the wrapped normal and von Mises processes are included, which can be seen as toroidal analogues of the Ornstein-Uhlenbeck diffusion. A collection of methods for approximate maximum likelihood estimation, organized in four blocks, is given: (i) based on the exact transition probability density, obtained as the numerical solution to the Fokker-Plank equation; (ii) based on wrapped pseudo-likelihoods; (iii) based on specific analytic approximations by wrapped processes; (iv) based on maximum likelihood of the stationary densities. The package allows the replicability of the results in García-Portugués et al. (2019) <doi:10.1007/s11222-017-9790-2>. |
| Authors: | Eduardo García-Portugués [aut, cre]
|
| Maintainer: | Eduardo García-Portugués <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.10 |
| Built: | 2024-10-30 06:52:23 UTC |
| Source: | CRAN |
Constructs drift matrices such that
solve(A) %*% Sigma is symmetric.
alphaToA(alpha, sigma = NULL, rho = 0, Sigma = NULL) aToAlpha(A, sigma = NULL, rho = 0, Sigma = NULL)alphaToA(alpha, sigma = NULL, rho = 0, Sigma = NULL) aToAlpha(A, sigma = NULL, rho = 0, Sigma = NULL)
alpha |
vector of length |
sigma |
vector of length |
rho |
correlation of |
Sigma |
the diffusion matrix of size |
A |
matrix of size |
The parametrization enforces that solve(A) %*% Sigma
is symmetric. Positive definiteness is guaranteed if alpha[3]^2 <
rho^2 * (alpha[1] - alpha[2])^2 / 4 + alpha[1] * alpha[2].
The drift matrix A or the alpha vector.
# Parameters alpha <- 3:1 Sigma <- rbind(c(1, 0.5), c(0.5, 4)) # Covariance matrix A <- alphaToA(alpha = alpha, Sigma = Sigma) S <- 0.5 * solve(A) %*% Sigma det(S) # Check aToAlpha(A = alphaToA(alpha = alpha, Sigma = Sigma), Sigma = Sigma) alphaToA(alpha = aToAlpha(A = A, Sigma = Sigma), Sigma = Sigma)# Parameters alpha <- 3:1 Sigma <- rbind(c(1, 0.5), c(0.5, 4)) # Covariance matrix A <- alphaToA(alpha = alpha, Sigma = Sigma) S <- 0.5 * solve(A) %*% Sigma det(S) # Check aToAlpha(A = alphaToA(alpha = alpha, Sigma = Sigma), Sigma = Sigma) alphaToA(alpha = aToAlpha(A = A, Sigma = Sigma), Sigma = Sigma)
Approximate Maximum Likelihood Estimation (MLE) for the Wrapped Normal (WN) in 1D using the wrapped Ornstein–Uhlenbeck diffusion.
approxMleWn1D(data, delta, start, alpha = NA, mu = NA, sigma = NA, lower = c(0.01, -pi, 0.01), upper = c(25, pi, 25), vmApprox = FALSE, maxK = 2, ...)approxMleWn1D(data, delta, start, alpha = NA, mu = NA, sigma = NA, lower = c(0.01, -pi, 0.01), upper = c(25, pi, 25), vmApprox = FALSE, maxK = 2, ...)
data |
a matrix of dimension |
delta |
discretization step. |
start |
starting values, a matrix with |
alpha, mu, sigma
|
if their values are provided, the likelihood function
is optimized with respect to the rest of unspecified parameters. The number
of elements in |
lower, upper
|
bound for box constraints as in method |
vmApprox |
flag to indicate von Mises approximation to wrapped normal.
See |
maxK |
maximum absolute winding number used if |
... |
further parameters passed to |
See Section 3.3 in García-Portugués et al. (2019) for details.
Output from mleOptimWrapper.
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
alpha <- 0.5 mu <- 0 sigma <- 2 samp <- rTrajWn1D(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 1000, delta = 0.1) approxMleWn1D(data = samp, delta = 0.1, start = c(alpha, mu, sigma)) approxMleWn1D(data = samp, delta = 0.1, sigma = sigma, start = c(alpha, mu), lower = c(0.01, -pi), upper = c(25, pi)) approxMleWn1D(data = samp, delta = 0.1, mu = mu, start = c(alpha, sigma), lower = c(0.01, 0.01), upper = c(25, 25))alpha <- 0.5 mu <- 0 sigma <- 2 samp <- rTrajWn1D(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 1000, delta = 0.1) approxMleWn1D(data = samp, delta = 0.1, start = c(alpha, mu, sigma)) approxMleWn1D(data = samp, delta = 0.1, sigma = sigma, start = c(alpha, mu), lower = c(0.01, -pi), upper = c(25, pi)) approxMleWn1D(data = samp, delta = 0.1, mu = mu, start = c(alpha, sigma), lower = c(0.01, 0.01), upper = c(25, 25))
Approximate Maximum Likelihood Estimation (MLE) for the Wrapped Normal (WN) in 2D using the wrapped Ornstein–Uhlenbeck diffusion.
approxMleWn2D(data, delta, start, alpha = rep(NA, 3), mu = rep(NA, 2), sigma = rep(NA, 2), rho = NA, lower = c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01, -0.99), upper = c(rep(25, 3), pi, pi, 25, 25, 0.99), maxK = 2, ...)approxMleWn2D(data, delta, start, alpha = rep(NA, 3), mu = rep(NA, 2), sigma = rep(NA, 2), rho = NA, lower = c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01, -0.99), upper = c(rep(25, 3), pi, pi, 25, 25, 0.99), maxK = 2, ...)
data |
a matrix of dimension |
delta |
discretization step. |
start |
starting values, a matrix with |
alpha, mu, sigma, rho
|
if their values are provided, the likelihood
function is optimized with respect to the rest of unspecified parameters.
The number of elements in |
lower, upper
|
bound for box constraints as in method |
maxK |
maximum absolute winding number used if |
... |
further parameters passed to |
See Section 3.3 in García-Portugués et al. (2019) for details.
Output from mleOptimWrapper.
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
alpha <- c(2, 2, -0.5) mu <- c(0, 0) sigma <- c(1, 1) rho <- 0.2 samp <- rTrajWn2D(x0 = c(0, 0), alpha = alpha, mu = mu, sigma = sigma, rho = rho, N = 1000, delta = 0.1) approxMleWn2D(data = samp, delta = 0.1, start = c(alpha, mu, sigma, rho)) approxMleWn2D(data = samp, delta = 0.1, alpha = alpha, start = c(mu, sigma), lower = c(-pi, -pi, 0.01, 0.01), upper = c(pi, pi, 25, 25)) mleMou(data = samp, delta = 0.1, start = c(alpha, mu, sigma), optMethod = "Nelder-Mead")alpha <- c(2, 2, -0.5) mu <- c(0, 0) sigma <- c(1, 1) rho <- 0.2 samp <- rTrajWn2D(x0 = c(0, 0), alpha = alpha, mu = mu, sigma = sigma, rho = rho, N = 1000, delta = 0.1) approxMleWn2D(data = samp, delta = 0.1, start = c(alpha, mu, sigma, rho)) approxMleWn2D(data = samp, delta = 0.1, alpha = alpha, start = c(mu, sigma), lower = c(-pi, -pi, 0.01, 0.01), upper = c(pi, pi, 25, 25)) mleMou(data = samp, delta = 0.1, start = c(alpha, mu, sigma), optMethod = "Nelder-Mead")
Approximate Maximum Likelihood Estimation (MLE) for the Wrapped Normal (WN) diffusion, using the wrapped Ornstein–Uhlenbeck diffusion and assuming initial stationarity.
approxMleWnPairs(data, delta, start = c(0, 0, 1, 1, 0, 1, 1), alpha = rep(NA, 3), mu = rep(NA, 2), sigma = rep(NA, 2), rho = NA, lower = c(-pi, -pi, 0.01, 0.01, -25, 0.01, 0.01, -0.99), upper = c(pi, pi, 25, 25, 25, 25, 25, 0.99), maxK = 2, expTrc = 30, ...)approxMleWnPairs(data, delta, start = c(0, 0, 1, 1, 0, 1, 1), alpha = rep(NA, 3), mu = rep(NA, 2), sigma = rep(NA, 2), rho = NA, lower = c(-pi, -pi, 0.01, 0.01, -25, 0.01, 0.01, -0.99), upper = c(pi, pi, 25, 25, 25, 25, 25, 0.99), maxK = 2, expTrc = 30, ...)
data |
a matrix of dimension |
delta |
discretization step. |
start |
starting values, a matrix with |
alpha |
vector of length |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
lower, upper
|
bound for box constraints as in method |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
... |
further parameters passed to |
Output from mleOptimWrapper.
mu <- c(0, 0) alpha <- c(1, 2, 0.5) sigma <- c(1, 1) rho <- 0.5 set.seed(4567345) begin <- rStatWn2D(n = 200, mu = mu, alpha = alpha, sigma = sigma) end <- t(apply(begin, 1, function(x) rTrajWn2D(x0 = x, alpha = alpha, mu = mu, sigma = sigma, rho = rho, N = 1, delta = 0.1)[2, ])) data <- cbind(begin, end) approxMleWnPairs(data = data, delta = 0.1, start = c(2, pi/2, 2, 0.5, 0, 2, 1, 0.5))mu <- c(0, 0) alpha <- c(1, 2, 0.5) sigma <- c(1, 1) rho <- 0.5 set.seed(4567345) begin <- rStatWn2D(n = 200, mu = mu, alpha = alpha, sigma = sigma) end <- t(apply(begin, 1, function(x) rTrajWn2D(x0 = x, alpha = alpha, mu = mu, sigma = sigma, rho = rho, N = 1, delta = 0.1)[2, ])) data <- cbind(begin, end) approxMleWnPairs(data = data, delta = 0.1, start = c(2, pi/2, 2, 0.5, 0, 2, 1, 0.5))
Implementation of the Crank–Nicolson scheme for solving the Fokker–Planck equation
where is the transition probability density of the circular diffusion
.
crankNicolson1D(u0, b, sigma2, N, deltat, Mx, deltax, imposePositive = 0L)crankNicolson1D(u0, b, sigma2, N, deltat, Mx, deltax, imposePositive = 0L)
u0 |
matrix of size |
b |
vector of length |
sigma2 |
vector of length |
N |
increasing integer vector of length |
deltat |
time step. |
Mx |
size of the equispaced spatial grid in |
deltax |
space grid discretization. |
imposePositive |
flag to indicate whether the solution should be transformed in order to be always larger than a given tolerance. This prevents spurious negative values. The tolerance will be taken as |
The function makes use of solvePeriodicTridiag for obtaining implicitly the next step in time of the solution.
If imposePositive = TRUE, the code implicitly assumes that the solution integrates to one at any step. This might b unrealistic if the initial condition is not properly represented in the grid (for example, highly concentrated density in a sparse grid).
If nt > 1, a matrix of size c(Mx, nt) containing the discretized solution at the required times.
If nt == 1, a matrix of size c(Mx, nu0) containing the discretized solution at a fixed time for different starting values.
Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Springer, New York. doi:10.1007/978-1-4899-7278-1
# Parameters Mx <- 200 N <- 200 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] times <- seq(0, 1, l = N + 1) u0 <- dWn1D(x, pi/2, 0.05) b <- driftWn1D(x, alpha = 1, mu = pi, sigma = 1) sigma2 <- rep(1, Mx) # Full trajectory of the solution (including initial condition) u <- crankNicolson1D(u0 = cbind(u0), b = b, sigma2 = sigma2, N = 0:N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx) # Mass conservation colMeans(u) * 2 * pi # Visualization of tpd plotSurface2D(times, x, z = t(u), levels = seq(0, 3, l = 50)) # Only final time v <- crankNicolson1D(u0 = cbind(u0), b = b, sigma2 = sigma2, N = N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx) sum(abs(u[, N + 1] - v))# Parameters Mx <- 200 N <- 200 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] times <- seq(0, 1, l = N + 1) u0 <- dWn1D(x, pi/2, 0.05) b <- driftWn1D(x, alpha = 1, mu = pi, sigma = 1) sigma2 <- rep(1, Mx) # Full trajectory of the solution (including initial condition) u <- crankNicolson1D(u0 = cbind(u0), b = b, sigma2 = sigma2, N = 0:N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx) # Mass conservation colMeans(u) * 2 * pi # Visualization of tpd plotSurface2D(times, x, z = t(u), levels = seq(0, 3, l = 50)) # Only final time v <- crankNicolson1D(u0 = cbind(u0), b = b, sigma2 = sigma2, N = N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx) sum(abs(u[, N + 1] - v))
Implementation of the Crank–Nicolson scheme for solving the Fokker–Planck equation
where is the transition probability density of the toroidal diffusion
crankNicolson2D(u0, bx, by, sigma2x, sigma2y, sigmaxy, N, deltat, Mx, deltax, My, deltay, imposePositive = 0L)crankNicolson2D(u0, bx, by, sigma2x, sigma2y, sigmaxy, N, deltat, Mx, deltax, My, deltay, imposePositive = 0L)
u0 |
matrix of size |
bx, by
|
matrices of size |
sigma2x, sigma2y, sigmaxy
|
matrices of size |
N |
increasing integer vector of length |
deltat |
time step. |
Mx, My
|
sizes of the equispaced spatial grids in |
deltax, deltay
|
space grid discretizations for each component. |
imposePositive |
flag to indicate whether the solution should be transformed in order to be always larger than a given tolerance. This prevents spurious negative values. The tolerance will be taken as |
The function makes use of solvePeriodicTridiag for obtaining implicitly the next step in time of the solution.
If imposePositive = TRUE, the code implicitly assumes that the solution integrates to one at any step. This might b unrealistic if the initial condition is not properly represented in the grid (for example, highly concentrated density in a sparse grid).
If nt > 1, a matrix of size c(Mx * My, nt) containing the discretized solution at the required times with the c(Mx, My) matrix stored column-wise.
If nt == 1, a matrix of size c(Mx * My, nu0) containing the discretized solution at a fixed time for different starting values.
Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Springer, New York. doi:10.1007/978-1-4899-7278-1
# Parameters Mx <- 100 My <- 100 N <- 200 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] y <- seq(-pi, pi, l = My + 1)[-c(My + 1)] m <- c(pi / 2, pi) p <- c(0, 1) u0 <- c(outer(dWn1D(x, p[1], 0.5), dWn1D(y, p[2], 0.5))) bx <- outer(x, y, function(x, y) 5 * sin(m[1] - x)) by <- outer(x, y, function(x, y) 5 * sin(m[2] - y)) sigma2 <- matrix(1, nrow = Mx, ncol = My) sigmaxy <- matrix(0.5, nrow = Mx, ncol = My) # Full trajectory of the solution (including initial condition) u <- crankNicolson2D(u0 = cbind(u0), bx = bx, by = by, sigma2x = sigma2, sigma2y = sigma2, sigmaxy = sigmaxy, N = 0:N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx, My = My, deltay = 2 * pi / My) # Mass conservation colMeans(u) * 4 * pi^2 # Only final time v <- crankNicolson2D(u0 = cbind(u0), bx = bx, by = by, sigma2x = sigma2, sigma2y = sigma2, sigmaxy = sigmaxy, N = N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx, My = My, deltay = 2 * pi / My) sum(abs(u[, N + 1] - v)) ## Not run: # Visualization of tpd library(manipulate) manipulate({ plotSurface2D(x, y, z = matrix(u[, j + 1], Mx, My), main = round(mean(u[, j + 1]) * 4 * pi^2, 4), levels = seq(0, 2, l = 21)) points(p[1], p[2], pch = 16) points(m[1], m[2], pch = 16) }, j = slider(0, N)) ## End(Not run)# Parameters Mx <- 100 My <- 100 N <- 200 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] y <- seq(-pi, pi, l = My + 1)[-c(My + 1)] m <- c(pi / 2, pi) p <- c(0, 1) u0 <- c(outer(dWn1D(x, p[1], 0.5), dWn1D(y, p[2], 0.5))) bx <- outer(x, y, function(x, y) 5 * sin(m[1] - x)) by <- outer(x, y, function(x, y) 5 * sin(m[2] - y)) sigma2 <- matrix(1, nrow = Mx, ncol = My) sigmaxy <- matrix(0.5, nrow = Mx, ncol = My) # Full trajectory of the solution (including initial condition) u <- crankNicolson2D(u0 = cbind(u0), bx = bx, by = by, sigma2x = sigma2, sigma2y = sigma2, sigmaxy = sigmaxy, N = 0:N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx, My = My, deltay = 2 * pi / My) # Mass conservation colMeans(u) * 4 * pi^2 # Only final time v <- crankNicolson2D(u0 = cbind(u0), bx = bx, by = by, sigma2x = sigma2, sigma2y = sigma2, sigmaxy = sigmaxy, N = N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx, My = My, deltay = 2 * pi / My) sum(abs(u[, N + 1] - v)) ## Not run: # Visualization of tpd library(manipulate) manipulate({ plotSurface2D(x, y, z = matrix(u[, j + 1], Mx, My), main = round(mean(u[, j + 1]) * 4 * pi^2, 4), levels = seq(0, 2, l = 21)) points(p[1], p[2], pch = 16) points(m[1], m[2], pch = 16) }, j = slider(0, N)) ## End(Not run)
Evaluation of the bivariate Sine von Mises density and its normalizing constant.
dBvm(x, mu, kappa, logConst = NULL) constBvm(M = 25, kappa)dBvm(x, mu, kappa, logConst = NULL) constBvm(M = 25, kappa)
x |
a matrix of size |
mu |
two-dimensional vector of circular means. |
kappa |
three-dimensional vector with concentrations
|
logConst |
logarithm of the normalizing constant. Computed if
|
M |
number of terms considered in the series expansion used for evaluating the normalizing constant. |
If or and ,
then constBvm will perform a Monte Carlo integration of the constant.
A vector of length nx with the evaluated density
(dBvm) or a scalar with the normaalizing constant (constBvm).
Singh, H., Hnizdo, V. and Demchuk, E. (2002) Probabilistic model for two dependent circular variables, Biometrika, 89(3):719–723, doi:10.1093/biomet/89.3.719
x <- seq(-pi, pi, l = 101)[-101] plotSurface2D(x, x, f = function(x) dBvm(x = x, mu = c(0, pi / 2), kappa = c(2, 3, 1)), fVect = TRUE)x <- seq(-pi, pi, l = 101)[-101] plotSurface2D(x, x, f = function(x) dBvm(x = x, mu = c(0, pi / 2), kappa = c(2, 3, 1)), fVect = TRUE)
Returns suitably lagged and iterated circular differences.
diffCirc(x, circular = TRUE, ...)diffCirc(x, circular = TRUE, ...)
x |
wrapped or unwrapped angles to be differenced. Must be a vector or a matrix, see details. |
circular |
convenience flag to indicate whether wrapping should be
done. If |
... |
parameters to be passed to |
If x is a matrix then the difference operations are carried
out row-wise, on each column separately.
The value of diff(x, ...), circularly wrapped. Default
parameters give an object of the kind of x with one less entry or row.
# Vectors x <- c(-pi, -pi/2, pi - 0.1, -pi + 0.2) diffCirc(x) - diff(x) # Matrices set.seed(234567) N <- 100 x <- t(euler2D(x0 = rbind(c(0, 0)), A = diag(c(1, 1)), sigma = rep(2, 2), mu = c(pi, pi), N = N, delta = 1, type = 2)[1, , ]) diffCirc(x) - diff(x)# Vectors x <- c(-pi, -pi/2, pi - 0.1, -pi + 0.2) diffCirc(x) - diff(x) # Matrices set.seed(234567) N <- 100 x <- t(euler2D(x0 = rbind(c(0, 0)), A = diag(c(1, 1)), sigma = rep(2, 2), mu = c(pi, pi), N = N, delta = 1, type = 2)[1, , ]) diffCirc(x) - diff(x)
Computes the circular density of Jones and Pewsey (2005).
dJp(x, mu, kappa, psi, const = NULL) constJp(mu, kappa, psi, M = 200)dJp(x, mu, kappa, psi, const = NULL) constJp(mu, kappa, psi, M = 200)
x |
evaluation angles, not necessary in |
mu |
circular mean. |
kappa |
non-negative concentration parameter. |
psi |
shape parameter, see details. |
const |
normalizing constant, computed with |
M |
grid size for computing the normalizing constant by numerical integration. |
Particular interesting choices for the shape parameter are:
psi = -1: gives the Wrapped Cauchy as stationary density.
psi = 0: is the sinusoidal drift of the vM diffusion.
psi = 1: gives the Cardioid as stationary density.
A vector of the same length as x containing the density.
Jones, M. C. and Pewsey, A. (2005). A family of symmetric distributions on the circle. Journal of the American Statistical Association, 100(472):1422–1428. doi:10.1198/016214505000000286
x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "Density", ylim = c(0, 0.6)) for (i in 0:20) { lines(x, dJp(x = x, mu = 0, kappa = 1, psi = -2 + 4 * i / 20), col = rainbow(21)[i + 1]) }x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "Density", ylim = c(0, 0.6)) for (i in 0:20) { lines(x, dJp(x = x, mu = 0, kappa = 1, psi = -2 + 4 * i / 20), col = rainbow(21)[i + 1]) }
Wrapped pseudo-transition probability densities.
dPsTpd(x, x0, t, method = c("E", "SO", "SO2"), b, jac.b, sigma2, b1, b2, circular = TRUE, maxK = 2, vmApprox = FALSE, twokpi = NULL, ...)dPsTpd(x, x0, t, method = c("E", "SO", "SO2"), b, jac.b, sigma2, b1, b2, circular = TRUE, maxK = 2, vmApprox = FALSE, twokpi = NULL, ...)
x |
a matrix of dimension |
x0 |
a matrix of dimension |
t |
time step between |
method |
a string for choosing |
b |
drift function. Must return a matrix of the same size as |
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 |
b1 |
first derivative of the drift function (univariate). Must return
a vector of the same length as |
b2 |
second derivative of the drift function (univariate). Must return
a vector of the same length as |
circular |
flag to indicate circular data. |
maxK |
maximum absolute winding number used if |
vmApprox |
flag to indicate von Mises approximation to wrapped normal.
See |
twokpi |
optional matrix of winding numbers to avoid its recomputation. See details. |
... |
additional parameters passed to |
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
andas.matrix(do.call(what = expand.grid,
args = rep(list(2 * pi * c(-maxK:maxK)), p))) otherwise.
Output from mleOptimWrapper.
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
# 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))# 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))
Drift for the Langevin diffusion associated to the Jones and Pewsey (JP) family of circular distributions.
driftJp(x, alpha, mu, psi)driftJp(x, alpha, mu, psi)
x |
vector with the evaluation points for the drift. |
alpha |
strength of the drift. |
mu |
unconditional mean of the diffusion. |
psi |
shape parameter, see details. |
Particular interesting choices for the shape parameter are:
psi = -1: gives the Wrapped Cauchy as stationary density.
psi = 0: is the sinusoidal drift of the vM diffusion.
psi = 1: gives the Cardioid as stationary density.
See Section 2.2.3 in García-Portugués et al. (2019) for details.
A vector of the same length as x containing the drift.
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
Jones, M. C. and Pewsey, A. (2005). A family of symmetric distributions on the circle. Journal of the American Statistical Association, 100(472):1422–1428. doi:10.1198/016214505000000286
x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 0:20) { lines(x, driftJp(x = x, alpha = 1, mu = 0, psi = -1 + 2 * i / 20), col = rainbow(21)[i + 1]) }x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 0:20) { lines(x, driftJp(x = x, alpha = 1, mu = 0, psi = -1 + 2 * i / 20), col = rainbow(21)[i + 1]) }
Drift for the Langevin diffusion associated to a mixture of
m independent (multivariate) von Mises (mivM) of dimension p.
driftMixIndVm(x, A, M, sigma, p, expTrc = 30)driftMixIndVm(x, A, M, sigma, p, expTrc = 30)
x |
matrix of size |
A |
matrix of size |
M |
matrix of size |
sigma |
diffusion coefficient. |
p |
vector of length |
expTrc |
truncation for exponential: |
driftMixVm is more efficient for the circular case.
The diffusion matrix is . See Section 2.2.4 in
García-Portugués et al. (2019) for details.
A matrix of the same size as x containing the drift.
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
# 1D x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 1:10) { lines(x, driftMixIndVm(x = cbind(x), A = cbind(c(2, 2)), M = cbind(c(0, -pi + 2 * pi * i / 10)), sigma = 1, p = c(0.5, 0.5)), col = rainbow(10)[i]) } # 2D x <- seq(-pi, pi, l = 100) plotSurface2D(x, x, f = function(x) sqrt(rowSums(driftMixIndVm(x = x, A = rbind(c(1, 1), c(1, 1)), M = rbind(c(1, 1), c(-1, -1)), sigma = 1, p = c(0.25, 0.75))^2)), fVect = TRUE)# 1D x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 1:10) { lines(x, driftMixIndVm(x = cbind(x), A = cbind(c(2, 2)), M = cbind(c(0, -pi + 2 * pi * i / 10)), sigma = 1, p = c(0.5, 0.5)), col = rainbow(10)[i]) } # 2D x <- seq(-pi, pi, l = 100) plotSurface2D(x, x, f = function(x) sqrt(rowSums(driftMixIndVm(x = x, A = rbind(c(1, 1), c(1, 1)), M = rbind(c(1, 1), c(-1, -1)), sigma = 1, p = c(0.25, 0.75))^2)), fVect = TRUE)
Drift for the Langevin diffusion associated to a mixture of
m independent von Mises (mivM) of dimension one.
driftMixVm(x, alpha, mu, sigma, p, expTrc = 30)driftMixVm(x, alpha, mu, sigma, p, expTrc = 30)
x |
vector with the evaluation points for the drift. |
alpha |
vector of length |
mu |
vector of length |
sigma |
diffusion coefficient. |
p |
vector of length |
expTrc |
truncation for exponential: |
driftMixIndVm is more general, but less efficient for
the circular case. See Section 2.2.4 in García-Portugués et al. (2019) for
details.
A vector of the same length as x containing the drift.
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
x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 1:10) { lines(x, driftMixVm(x = x, alpha = c(2, 2), mu = c(0, -pi + 2 * pi * i / 10), sigma = 1, p = c(0.5, 0.5)), col = rainbow(10)[i]) }x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 1:10) { lines(x, driftMixVm(x = x, alpha = c(2, 2), mu = c(0, -pi + 2 * pi * i / 10), sigma = 1, p = c(0.5, 0.5)), col = rainbow(10)[i]) }
Drift for the Langevin diffusion associated to the Multivariate
von Mises (MvM) in dimension p.
driftMvm(x, alpha, mu, A = 0)driftMvm(x, alpha, mu, A = 0)
x |
matrix of size |
alpha |
vector of length |
mu |
vector of length |
A |
matrix of size |
See Section 2.2.1 in García-Portugués et al. (2019) for details.
A matrix of the same size as x containing the drift.
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
# 1D x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 0:20) { lines(x, driftMvm(x = x, alpha = 3 * i / 20, mu = 0, A = 0), col = rainbow(21)[i + 1]) } # 2D x <- seq(-pi, pi, l = 100) plotSurface2D(x, x, f = function(x) sqrt(rowSums(driftMvm(x = x, alpha = c(2, 2), mu = c(-1, -1), A = rbind(c(0, 0), c(0, 0)))^2)), fVect = TRUE)# 1D x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 0:20) { lines(x, driftMvm(x = x, alpha = 3 * i / 20, mu = 0, A = 0), col = rainbow(21)[i + 1]) } # 2D x <- seq(-pi, pi, l = 100) plotSurface2D(x, x, f = function(x) sqrt(rowSums(driftMvm(x = x, alpha = c(2, 2), mu = c(-1, -1), A = rbind(c(0, 0), c(0, 0)))^2)), fVect = TRUE)
Drift for the Langevin diffusion associated to the
(multivariate) Wrapped Normal (WN) in dimension p.
driftWn(x, A, mu, Sigma, invSigmaA = NULL, maxK = 2, expTrc = 30)driftWn(x, A, mu, Sigma, invSigmaA = NULL, maxK = 2, expTrc = 30)
x |
matrix of size |
A |
matrix of size |
mu |
vector of length |
Sigma |
diffusion matrix, of size |
invSigmaA |
the matrix |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
See Section 2.2.2 in García-Portugués et al. (2019) for details.
driftWn1D and driftWn2D are more
efficient for the 1D and 2D cases.
A matrix of the same size as x containing the drift.
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
# 1D x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 1:20) { lines(x, driftWn(x = cbind(x), A = 1 * i / 20, mu = 0, Sigma = 1), col = rainbow(20)[i]) } # 2D x <- seq(-pi, pi, l = 100) plotSurface2D(x, x, f = function(x) sqrt(rowSums( driftWn(x = x, A = alphaToA(alpha = c(1, 1, 0.5), sigma = c(1.5, 1.5)), mu = c(0, 0), Sigma = diag(c(1.5^2, 1.5^2)))^2)), fVect = TRUE)# 1D x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "drift") for (i in 1:20) { lines(x, driftWn(x = cbind(x), A = 1 * i / 20, mu = 0, Sigma = 1), col = rainbow(20)[i]) } # 2D x <- seq(-pi, pi, l = 100) plotSurface2D(x, x, f = function(x) sqrt(rowSums( driftWn(x = x, A = alphaToA(alpha = c(1, 1, 0.5), sigma = c(1.5, 1.5)), mu = c(0, 0), Sigma = diag(c(1.5^2, 1.5^2)))^2)), fVect = TRUE)
Computes the drift of the WN diffusion in 1D in a vectorized way.
driftWn1D(x, alpha, mu, sigma, maxK = 2L, expTrc = 30)driftWn1D(x, alpha, mu, sigma, maxK = 2L, expTrc = 30)
x |
a vector of length |
alpha |
drift parameter. |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
A vector of length n containing the drift evaluated at x.
driftWn1D(x = seq(0, pi, l = 10), alpha = 1, mu = 0, sigma = 1, maxK = 2, expTrc = 30)driftWn1D(x = seq(0, pi, l = 10), alpha = 1, mu = 0, sigma = 1, maxK = 2, expTrc = 30)
Computes the drift of the WN diffusion in 2D in a vectorized way.
driftWn2D(x, A, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)driftWn2D(x, A, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)
x |
a matrix of dimension |
A |
drift matrix of size |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
A matrix of size c(n, 2) containing the drift evaluated at x.
alpha <- 3:1 mu <- c(0, 0) sigma <- 1:2 rho <- 0.5 Sigma <- diag(sigma^2) Sigma[1, 2] <- Sigma[2, 1] <- rho * prod(sigma) A <- alphaToA(alpha = alpha, sigma = sigma, rho = rho) x <- rbind(c(0, 1), c(1, 0.1), c(pi, pi), c(-pi, -pi), c(pi / 2, 0)) driftWn2D(x = x, A = A, mu = mu, sigma = sigma, rho = rho) driftWn(x = x, A = A, mu = c(0, 0), Sigma = Sigma)alpha <- 3:1 mu <- c(0, 0) sigma <- 1:2 rho <- 0.5 Sigma <- diag(sigma^2) Sigma[1, 2] <- Sigma[2, 1] <- rho * prod(sigma) A <- alphaToA(alpha = alpha, sigma = sigma, rho = rho) x <- rbind(c(0, 1), c(1, 0.1), c(pi, pi), c(-pi, -pi), c(pi / 2, 0)) driftWn2D(x = x, A = A, mu = mu, sigma = sigma, rho = rho) driftWn(x = x, A = A, mu = c(0, 0), Sigma = Sigma)
Stationary density of the WN diffusion.
dStatWn2D(x, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)dStatWn2D(x, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)
x |
a matrix of dimension |
alpha |
vector of length |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
A vector of size n containing the stationary density evaluated at x.
set.seed(345567) alpha <- c(2, 1, -1) sigma <- c(1.5, 2) Sigma <- diag(sigma^2) A <- alphaToA(alpha = alpha, sigma = sigma) mu <- c(pi, pi) dStatWn2D(x = toPiInt(matrix(1:20, nrow = 10, ncol = 2)), mu = mu, alpha = alpha, sigma = sigma) dTpdWou(t = 10, x = toPiInt(matrix(1:20, nrow = 10, ncol = 2)), A = A, mu = mu, Sigma = Sigma, x0 = mu) xth <- seq(-pi, pi, l = 100) contour(xth, xth, matrix(dStatWn2D(x = as.matrix(expand.grid(xth, xth)), alpha = alpha, sigma = sigma, mu = mu), nrow = length(xth), ncol = length(xth)), nlevels = 50) points(rStatWn2D(n = 1000, mu = mu, alpha = alpha, sigma = sigma), col = 2)set.seed(345567) alpha <- c(2, 1, -1) sigma <- c(1.5, 2) Sigma <- diag(sigma^2) A <- alphaToA(alpha = alpha, sigma = sigma) mu <- c(pi, pi) dStatWn2D(x = toPiInt(matrix(1:20, nrow = 10, ncol = 2)), mu = mu, alpha = alpha, sigma = sigma) dTpdWou(t = 10, x = toPiInt(matrix(1:20, nrow = 10, ncol = 2)), A = A, mu = mu, Sigma = Sigma, x0 = mu) xth <- seq(-pi, pi, l = 100) contour(xth, xth, matrix(dStatWn2D(x = as.matrix(expand.grid(xth, xth)), alpha = alpha, sigma = sigma, mu = mu), nrow = length(xth), ncol = length(xth)), nlevels = 50) points(rStatWn2D(n = 1000, mu = mu, alpha = alpha, sigma = sigma), col = 2)
Transition probability density of the multivariate Ornstein–Uhlenbeck (OU) diffusion
dTpdMou(x, x0, t, A, mu, Sigma, eigA = NULL, log = FALSE) meantMou(t, x0, A, mu, eigA = NULL) covtMou(t, A, Sigma, eigA = NULL)dTpdMou(x, x0, t, A, mu, Sigma, eigA = NULL, log = FALSE) meantMou(t, x0, A, mu, eigA = NULL) covtMou(t, A, Sigma, eigA = NULL)
x |
matrix of with |
x0 |
initial point. |
t |
time between observations. |
A |
the drift matrix, of size |
mu |
unconditional mean of the diffusion. |
Sigma |
square of the diffusion matrix, a matrix of size |
eigA |
optional argument containing |
log |
flag to indicate whether to compute the logarithm of the density. |
The transition probability density is a multivariate normal with
mean meantMou and covariance covtMou. See
dTpdOu for the univariate case (more efficient).
solve(A) %*% Sigma has to be a covariance matrix (symmetric and
positive definite) in order to have a proper transition density. For the
bivariate case, this can be ensured with the alphaToA function.
In the multivariate case, it is ensured if Sigma is isotropic and
A is a covariance matrix.
A matrix of the same size as x containing the evaluation of
the density.
x <- seq(-4, 4, by = 0.1) xx <- as.matrix(expand.grid(x, x)) isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate( image(x, x, matrix(dTpdMou(x = xx, x0 = c(1, 2), t = t, A = alphaToA(alpha = c(1, 2, 0.5), sigma = 1:2), mu = c(0, 0), Sigma = diag((1:2)^2)), nrow = length(x), ncol = length(x)), zlim = c(0, 0.25)), t = manipulate::slider(0.1, 5, step = 0.1)) }x <- seq(-4, 4, by = 0.1) xx <- as.matrix(expand.grid(x, x)) isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate( image(x, x, matrix(dTpdMou(x = xx, x0 = c(1, 2), t = t, A = alphaToA(alpha = c(1, 2, 0.5), sigma = 1:2), mu = c(0, 0), Sigma = diag((1:2)^2)), nrow = length(x), ncol = length(x)), zlim = c(0, 0.25)), t = manipulate::slider(0.1, 5, step = 0.1)) }
Transition probability density of the univariate Ornstein–Uhlenbeck (OU) diffusion
dTpdOu(x, x0, t, alpha, mu, sigma, log = FALSE) meantOu(x0, t, alpha, mu) vartOu(t, alpha, sigma) covstOu(s, t, alpha, sigma)dTpdOu(x, x0, t, alpha, mu, sigma, log = FALSE) meantOu(x0, t, alpha, mu) vartOu(t, alpha, sigma) covstOu(s, t, alpha, sigma)
x |
vector with the evaluation points. |
x0 |
initial point. |
t, s
|
time between observations. |
alpha |
strength of the drift. |
mu |
unconditional mean of the diffusion. |
sigma |
diffusion coefficient. |
log |
flag to indicate whether to compute the logarithm of the density. |
The transition probability density is a normal density with mean
meantOu and variance vartOu. See
dTpdMou for the multivariate case (less efficient for
dimension one).
A vector of the same length as x containing the evaluation of
the density.
x <- seq(-4, 4, by = 0.01) plot(x, dTpdOu(x = x, x0 = 3, t = 0.1, alpha = 1, mu = -1, sigma = 1), type = "l", ylim = c(0, 1.5), xlab = "x", ylab = "Density", col = rainbow(20)[1]) for (i in 2:20) { lines(x, dTpdOu(x = x, x0 = 3, t = i / 10, alpha = 1, mu = -1, sigma = 1), col = rainbow(20)[i]) }x <- seq(-4, 4, by = 0.01) plot(x, dTpdOu(x = x, x0 = 3, t = 0.1, alpha = 1, mu = -1, sigma = 1), type = "l", ylim = c(0, 1.5), xlab = "x", ylab = "Density", col = rainbow(20)[1]) for (i in 2:20) { lines(x, dTpdOu(x = x, x0 = 3, t = i / 10, alpha = 1, mu = -1, sigma = 1), col = rainbow(20)[i]) }
Computation of the transition probability density (tpd) of the Wrapped Normal (WN) or von Mises (vM) diffusion, by solving its associated Fokker–Planck Partial Differential Equation (PDE) in 1D.
dTpdPde1D(Mx = 500, x0, t, alpha, mu, sigma, type = "WN", Mt = ceiling(100 * t), sdInitial = 0.1, ...)dTpdPde1D(Mx = 500, x0, t, alpha, mu, sigma, type = "WN", Mt = ceiling(100 * t), sdInitial = 0.1, ...)
Mx |
size of the equispaced spatial grid in |
x0 |
point giving the mean of the initial circular density, a WN with
standard deviation equal to |
t |
time separating |
alpha |
drift parameter. |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
type |
either |
Mt |
size of the time grid in |
sdInitial |
the standard deviation of the concentrated WN giving the initial condition. |
... |
Further parameters passed to |
A combination of small sdInitial and coarse space-time
discretization (small Mx and Mt) is prone to create numerical
instabilities. See Sections 3.4.1, 2.2.1 and 2.2.2 in García-Portugués et
al. (2019) for details.
A vector of length Mx with the tpd evaluated at
seq(-pi, pi, l = Mx + 1)[-(Mx + 1)].
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
Mx <- 100 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] x0 <- pi t <- 0.5 alpha <- 1 mu <- 0 sigma <- 1 isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate({ plot(x, dTpdPde1D(Mx = Mx, x0 = x0, t = t, alpha = alpha, mu = 0, sigma = sigma), type = "l", ylab = "Density", xlab = "", ylim = c(0, 0.75)) lines(x, dTpdWou1D(x = x, x0 = rep(x0, Mx), t = t, alpha = alpha, mu = 0, sigma = sigma), col = 2) }, x0 = manipulate::slider(-pi, pi, step = 0.01, initial = 0), alpha = manipulate::slider(0.01, 5, step = 0.01, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.01, initial = 1), t = manipulate::slider(0.01, 5, step = 0.01, initial = 1)) }Mx <- 100 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] x0 <- pi t <- 0.5 alpha <- 1 mu <- 0 sigma <- 1 isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate({ plot(x, dTpdPde1D(Mx = Mx, x0 = x0, t = t, alpha = alpha, mu = 0, sigma = sigma), type = "l", ylab = "Density", xlab = "", ylim = c(0, 0.75)) lines(x, dTpdWou1D(x = x, x0 = rep(x0, Mx), t = t, alpha = alpha, mu = 0, sigma = sigma), col = 2) }, x0 = manipulate::slider(-pi, pi, step = 0.01, initial = 0), alpha = manipulate::slider(0.01, 5, step = 0.01, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.01, initial = 1), t = manipulate::slider(0.01, 5, step = 0.01, initial = 1)) }
Computation of the transition probability density (tpd) of the Wrapped Normal (WN) or Multivariate von Mises (MvM) diffusion, by solving its associated Fokker–Planck Partial Differential Equation (PDE) in 2D.
dTpdPde2D(Mx = 50, My = 50, x0, t, alpha, mu, sigma, rho = 0, type = "WN", Mt = ceiling(100 * t), sdInitial = 0.1, ...)dTpdPde2D(Mx = 50, My = 50, x0, t, alpha, mu, sigma, rho = 0, type = "WN", Mt = ceiling(100 * t), sdInitial = 0.1, ...)
Mx, My
|
sizes of the equispaced spatial grids in |
x0 |
point giving the mean of the initial circular density, an
isotropic WN with standard deviations equal to |
t |
time separating |
alpha |
for |
mu |
vector of length |
sigma |
for |
rho |
for |
type |
either |
Mt |
size of the time grid in |
sdInitial |
standard deviations of the concentrated WN giving the initial condition. |
... |
Further parameters passed to |
A combination of small sdInitial and coarse space-time
discretization (small Mx and Mt) is prone to create numerical
instabilities. See Sections 3.4.2, 2.2.1 and 2.2.2 in García-Portugués et al.
(2019) for details.
A matrix of size c(Mx, My) with the tpd evaluated at the
combinations of seq(-pi, pi, l = Mx + 1)[-(Mx + 1)] and
seq(-pi, pi, l = My + 1)[-(My + 1)].
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
M <- 100 x <- seq(-pi, pi, l = M + 1)[-c(M + 1)] image(x, x, dTpdPde2D(Mx = M, My = M, x0 = c(0, pi), t = 1, alpha = c(1, 1, 0.5), mu = c(pi / 2, 0), sigma = 1:2), zlim = c(0, 0.25), col = matlab.like.colorRamps(20), xlab = "x", ylab = "y")M <- 100 x <- seq(-pi, pi, l = M + 1)[-c(M + 1)] image(x, x, dTpdPde2D(Mx = M, My = M, x0 = c(0, pi), t = 1, alpha = c(1, 1, 0.5), mu = c(pi / 2, 0), sigma = 1:2), zlim = c(0, 0.25), col = matlab.like.colorRamps(20), xlab = "x", ylab = "y")
Conditional probability density of the Wrapped Ornstein–Uhlenbeck (WOU) process.
dTpdWou(x, t, A, mu, Sigma, x0, maxK = 2, eigA = NULL, invASigma = NULL)dTpdWou(x, t, A, mu, Sigma, x0, maxK = 2, eigA = NULL, invASigma = NULL)
x |
matrix of size |
t |
a scalar containing the times separating |
A |
matrix of size |
mu |
mean parameter. Must be in |
Sigma |
diffusion matrix, of size |
x0 |
vector of length |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
eigA |
optional argument containing |
invASigma |
the matrix |
See Section 3.3 in García-Portugués et al. (2019) for details.
dTpdWou1D and dTpdWou2D are more efficient
implementations for the 1D and 2D cases, respectively.
A vector of length n with the density evaluated at x.
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
# 1D t <- 0.5 alpha <- 1 mu <- 0 sigma <- 1 x0 <- pi x <- seq(-pi, pi, l = 10) dTpdWou(x = cbind(x), x0 = x0, t = t, A = alpha, mu = 0, Sigma = sigma^2) - dTpdWou1D(x = cbind(x), x0 = rep(x0, 10), t = t, alpha = alpha, mu = 0, sigma = sigma) # 2D t <- 0.5 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) mu <- c(pi, 0) x0 <- c(0, 0) x <- seq(-pi, pi, l = 5) x <- as.matrix(expand.grid(x, x)) dTpdWou(x = x, x0 = x0, t = t, A = A, mu = mu, Sigma = Sigma) - dTpdWou2D(x = x, x0 = rbind(x0), t = t, alpha = alpha, mu = mu, sigma = sigma, rho = rho)# 1D t <- 0.5 alpha <- 1 mu <- 0 sigma <- 1 x0 <- pi x <- seq(-pi, pi, l = 10) dTpdWou(x = cbind(x), x0 = x0, t = t, A = alpha, mu = 0, Sigma = sigma^2) - dTpdWou1D(x = cbind(x), x0 = rep(x0, 10), t = t, alpha = alpha, mu = 0, sigma = sigma) # 2D t <- 0.5 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) mu <- c(pi, 0) x0 <- c(0, 0) x <- seq(-pi, pi, l = 5) x <- as.matrix(expand.grid(x, x)) dTpdWou(x = x, x0 = x0, t = t, A = A, mu = mu, Sigma = Sigma) - dTpdWou2D(x = x, x0 = rbind(x0), t = t, alpha = alpha, mu = mu, sigma = sigma, rho = rho)
Computation of the transition probability density (tpd) for a WN diffusion.
dTpdWou1D(x, x0, t, alpha, mu, sigma, maxK = 2L, expTrc = 30, vmApprox = 0L, kt = 0, logConstKt = 0)dTpdWou1D(x, x0, t, alpha, mu, sigma, maxK = 2L, expTrc = 30, vmApprox = 0L, kt = 0, logConstKt = 0)
x |
a vector of length |
x0 |
a vector of length |
t |
a scalar containing the times separating |
alpha |
drift parameter. |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
vmApprox |
whether to use the von Mises approximation to a wrapped normal ( |
kt |
concentration for the von Mises, a suitable output from |
logConstKt |
the logarithm of the von Mises normalizing constant associated to the concentration |
See Section 3.3 in García-Portugués et al. (2019) for details. See dTpdWou for the general case (less efficient for 2D).
A vector of size n containing the tpd evaluated at x.
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
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)))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)))
Computation of the transition probability density (tpd) for a WN diffusion (with diagonal diffusion matrix)
dTpdWou2D(x, x0, t, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)dTpdWou2D(x, x0, t, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)
x |
a matrix of dimension |
x0 |
a matrix of dimension |
t |
a scalar containing the times separating |
alpha |
vector of length |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
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).
A vector of size n containing the tpd evaluated at x.
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
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))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))
Computes the density of a von Mises in a numerically stable way.
dVm(x, mu, kappa)dVm(x, mu, kappa)
x |
evaluation angles, not necessary in |
mu |
circular mean. |
kappa |
non-negative concentration parameter. |
A vector of the same length as x containing the density.
Jammalamadaka, S. R. and SenGupta, A. (2001) Topics in Circular Statistics. World Scientific, Singapore. doi:10.1142/4031
x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "Density", ylim = c(0, 1)) for (i in 0:20) { lines(x, dVm(x = x, mu = 0, kappa = 5 * i / 20), col = rainbow(21)[i + 1]) }x <- seq(-pi, pi, l = 200) plot(x, x, type = "n", ylab = "Density", ylim = c(0, 1)) for (i in 0:20) { lines(x, dVm(x = x, mu = 0, kappa = 5 * i / 20), col = rainbow(21)[i + 1]) }
Computation of the WN density in 1D.
dWn1D(x, mu, sigma, maxK = 2L, expTrc = 30, vmApprox = 0L, kt = 0, logConstKt = 0)dWn1D(x, mu, sigma, maxK = 2L, expTrc = 30, vmApprox = 0L, kt = 0, logConstKt = 0)
x |
a vector of length |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
vmApprox |
whether to use the von Mises approximation to a wrapped normal ( |
kt |
concentration for the von Mises, a suitable output from |
logConstKt |
the logarithm of the von Mises normalizing constant associated to the concentration |
A vector of size n containing the density evaluated at x.
mu <- 0 sigma <- 1 dWn1D(x = seq(-pi, pi, l = 10), mu = mu, sigma = sigma, vmApprox = 0) # von Mises approximation kt <- scoreMatchWnVm(sigma2 = sigma^2) dWn1D(x = seq(-pi, pi, l = 10), mu = mu, sigma = sigma, vmApprox = 1, kt = kt, logConstKt = -log(2 * pi * besselI(x = kt, nu = 0, expon.scaled = TRUE)))mu <- 0 sigma <- 1 dWn1D(x = seq(-pi, pi, l = 10), mu = mu, sigma = sigma, vmApprox = 0) # von Mises approximation kt <- scoreMatchWnVm(sigma2 = sigma^2) dWn1D(x = seq(-pi, pi, l = 10), mu = mu, sigma = sigma, vmApprox = 1, kt = kt, logConstKt = -log(2 * pi * besselI(x = kt, nu = 0, expon.scaled = TRUE)))
Simulation of the Wrapped Normal (WN) diffusion or von Mises (vM) diffusion by the Euler method in 1D, for several starting values.
euler1D(x0, alpha, mu, sigma, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)euler1D(x0, alpha, mu, sigma, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)
x0 |
vector of length |
alpha |
drift parameter. |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
N |
number of discretization steps. |
delta |
discretization step. |
type |
integer giving the type of diffusion. Currently, only |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
A matrix of size c(nx0, N + 1) containing the nx0 discretized trajectories. The first column corresponds to the vector x0.
N <- 100 nx0 <- 20 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] set.seed(12345678) samp <- euler1D(x0 = x0, mu = 0, alpha = 3, sigma = 1, N = N, delta = 0.01, type = 2) tt <- seq(0, 1, l = N + 1) plot(rep(0, nx0), x0, pch = 16, col = rainbow(nx0), xlim = c(0, 1)) for (i in 1:nx0) linesCirc(tt, samp[i, ], col = rainbow(nx0)[i])N <- 100 nx0 <- 20 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] set.seed(12345678) samp <- euler1D(x0 = x0, mu = 0, alpha = 3, sigma = 1, N = N, delta = 0.01, type = 2) tt <- seq(0, 1, l = N + 1) plot(rep(0, nx0), x0, pch = 16, col = rainbow(nx0), xlim = c(0, 1)) for (i in 1:nx0) linesCirc(tt, samp[i, ], col = rainbow(nx0)[i])
Simulation of the Wrapped Normal (WN) diffusion or Multivariate von Mises (MvM) diffusion by the Euler method in 2D, for several starting values.
euler2D(x0, A, mu, sigma, rho = 0, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)euler2D(x0, A, mu, sigma, rho = 0, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)
x0 |
matrix of size |
A |
drift matrix of size |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
N |
number of discretization steps. |
delta |
discretization step. |
type |
integer giving the type of diffusion. Currently, only |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
An array of size c(nx0, 2, N + 1) containing the nx0 discretized trajectories. The first slice corresponds to the matrix x0.
N <- 100 nx0 <- 5 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] x0 <- as.matrix(expand.grid(x0, x0)) nx0 <- nx0^2 set.seed(12345678) samp <- euler2D(x0 = x0, mu = c(0, 0), A = rbind(c(3, 1), 1:2), sigma = c(1, 1), N = N, delta = 0.01, type = 2) plot(x0[, 1], x0[, 2], xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 16, col = rainbow(nx0)) for (i in 1:nx0) linesTorus(samp[i, 1, ], samp[i, 2, ], col = rainbow(nx0, alpha = 0.5)[i])N <- 100 nx0 <- 5 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] x0 <- as.matrix(expand.grid(x0, x0)) nx0 <- nx0^2 set.seed(12345678) samp <- euler2D(x0 = x0, mu = c(0, 0), A = rbind(c(3, 1), 1:2), sigma = c(1, 1), N = N, delta = 0.01, type = 2) plot(x0[, 1], x0[, 2], xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 16, col = rainbow(nx0)) for (i in 1:nx0) linesTorus(samp[i, 1, ], samp[i, 2, ], col = rainbow(nx0, alpha = 0.5)[i])
Joins the corresponding points with line segments or arrows that
exhibit wrapping in in the vertical axis.
linesCirc(x = seq_along(y), y, col = 1, lty = 1, ltyCross = lty, arrows = FALSE, ...)linesCirc(x = seq_along(y), y, col = 1, lty = 1, ltyCross = lty, arrows = FALSE, ...)
x |
vector with horizontal coordinates. |
y |
vector with vertical coordinates, wrapped in |
col |
color vector of length |
lty |
line type as in |
ltyCross |
specific line type for crossing segments. |
arrows |
flag for drawing arrows instead of line segments. |
... |
y is wrapped to before plotting.
Nothing. The functions are called for drawing wrapped lines.
x <- 1:100 y <- toPiInt(pi * cos(2 * pi * x / 100) + 0.5 * runif(50, -pi, pi)) plot(x, y, ylim = c(-pi, pi)) linesCirc(x = x, y = y, col = rainbow(length(x)), ltyCross = 2) plot(x, y, ylim = c(-pi, pi)) linesCirc(x = x, y = y, col = rainbow(length(x)), arrows = TRUE)x <- 1:100 y <- toPiInt(pi * cos(2 * pi * x / 100) + 0.5 * runif(50, -pi, pi)) plot(x, y, ylim = c(-pi, pi)) linesCirc(x = x, y = y, col = rainbow(length(x)), ltyCross = 2) plot(x, y, ylim = c(-pi, pi)) linesCirc(x = x, y = y, col = rainbow(length(x)), arrows = TRUE)
Joins the corresponding points with line segments or arrows that
exhibit wrapping in in the horizontal and vertical axes.
linesTorus(x, y, col = 1, lty = 1, ltyCross = lty, arrows = FALSE, ...)linesTorus(x, y, col = 1, lty = 1, ltyCross = lty, arrows = FALSE, ...)
x |
vector with horizontal coordinates, wrapped in |
y |
vector with vertical coordinates, wrapped in |
col |
color vector of length |
lty |
line type as in |
ltyCross |
specific line type for crossing segments. |
arrows |
flag for drawing arrows instead of line segments. |
... |
x and y are wrapped to before
plotting.
Nothing. The functions are called for drawing wrapped lines.
x <- toPiInt(rnorm(50, mean = seq(-pi, pi, l = 50), sd = 0.5)) y <- toPiInt(x + rnorm(50, mean = seq(-pi, pi, l = 50), sd = 0.5)) plot(x, y, xlim = c(-pi, pi), ylim = c(-pi, pi), col = rainbow(length(x)), pch = 19) linesTorus(x = x, y = y, col = rainbow(length(x)), ltyCross = 2) plot(x, y, xlim = c(-pi, pi), ylim = c(-pi, pi), col = rainbow(length(x)), pch = 19) linesTorus(x = x, y = y, col = rainbow(length(x)), arrows = TRUE)x <- toPiInt(rnorm(50, mean = seq(-pi, pi, l = 50), sd = 0.5)) y <- toPiInt(x + rnorm(50, mean = seq(-pi, pi, l = 50), sd = 0.5)) plot(x, y, xlim = c(-pi, pi), ylim = c(-pi, pi), col = rainbow(length(x)), pch = 19) linesTorus(x = x, y = y, col = rainbow(length(x)), ltyCross = 2) plot(x, y, xlim = c(-pi, pi), ylim = c(-pi, pi), col = rainbow(length(x)), pch = 19) linesTorus(x = x, y = y, col = rainbow(length(x)), arrows = TRUE)
Joins the corresponding points with line segments or arrows that
exhibit wrapping in in the horizontal and vertical axes.
linesTorus3d(x, y, z, col = 1, arrows = FALSE, ...)linesTorus3d(x, y, z, col = 1, arrows = FALSE, ...)
x, y
|
vectors with horizontal coordinates, wrapped in |
z |
vector with vertical coordinates, wrapped in |
col |
color vector of length |
arrows |
flag for drawing arrows instead of line segments. |
... |
x, y, and z are wrapped to
before plotting. arrows = TRUE makes sequential calls to
arrow3d, and is substantially slower than
arrows = FALSE.
Nothing. The functions are called for drawing wrapped lines.
if (requireNamespace("rgl")) { n <- 20 x <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) y <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) z <- toPiInt(x + y + rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) rgl::plot3d(x, y, z, xlim = c(-pi, pi), ylim = c(-pi, pi), zlim = c(-pi, pi), col = rainbow(n), size = 2, box = FALSE, axes = FALSE) linesTorus3d(x = x, y = y, z = z, col = rainbow(n), lwd = 2) torusAxis3d() rgl::plot3d(x, y, z, xlim = c(-pi, pi), ylim = c(-pi, pi), zlim = c(-pi, pi), col = rainbow(n), size = 2, box = FALSE, axes = FALSE) linesTorus3d(x = x, y = y, z = z, col = rainbow(n), ltyCross = 2, arrows = TRUE, theta = 0.1 * pi / 180, barblen = 0.1) torusAxis3d() }if (requireNamespace("rgl")) { n <- 20 x <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) y <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) z <- toPiInt(x + y + rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) rgl::plot3d(x, y, z, xlim = c(-pi, pi), ylim = c(-pi, pi), zlim = c(-pi, pi), col = rainbow(n), size = 2, box = FALSE, axes = FALSE) linesTorus3d(x = x, y = y, z = z, col = rainbow(n), lwd = 2) torusAxis3d() rgl::plot3d(x, y, z, xlim = c(-pi, pi), ylim = c(-pi, pi), zlim = c(-pi, pi), col = rainbow(n), size = 2, box = FALSE, axes = FALSE) linesTorus3d(x = x, y = y, z = z, col = rainbow(n), ltyCross = 2, arrows = TRUE, theta = 0.1 * pi / 180, barblen = 0.1) torusAxis3d() }
Computation of and the
inverse of .
logBesselI0Scaled(x, splineApprox = TRUE) a1Inv(x, splineApprox = TRUE)logBesselI0Scaled(x, splineApprox = TRUE) a1Inv(x, splineApprox = TRUE)
x |
evaluation vector. For |
splineApprox |
whether to use a pre-computed spline approximation (faster) or not. |
Both functions may rely on pre-computed spline interpolations
(logBesselI0ScaledSpline and a1InvSpline). Otherwise, a call
to besselI is done for and
is solved numerically. The data in which the
interpolation is based is given in the examples.
For x larger than 5e4, the asymptotic expansion of
besselIasym is employed.
A vector of the same length as x.
# Data employed for log besselI0 scaled x1 <- c(seq(0, 1, by = 1e-4), seq(1 + 1e-2, 10, by = 1e-3), seq(10 + 1e-1, 100, by = 1e-2), seq(100 + 1e0, 1e3, by = 1e0), seq(1000 + 1e1, 5e4, by = 2e1)) logBesselI0ScaledEvalGrid <- log(besselI(x = x1, nu = 0, expon.scaled = TRUE)) # save(list = "logBesselI0ScaledEvalGrid", # file = "logBesselI0ScaledEvalGrid.rda", compress = TRUE) # Data employed for A1 inverse x2 <- rev(c(seq(1e-04, 0.9 - 1e-4, by = 1e-4), seq(0.9, 1 - 1e-05, by = 1e-5))) a1InvEvalGrid <- sapply(x2, function(k) { uniroot(f = function(x) k - besselI(x, nu = 1, expon.scaled = TRUE) / besselI(x, nu = 0, expon.scaled = TRUE), lower = 1e-06, upper = 1e+05, tol = 1e-15)$root }) # save(list = "a1InvEvalGrid", file = "a1InvEvalGrid.rda", compress = TRUE) # Accuracy logBesselI0Scaled x <- seq(0, 1e3, l = 1e3) summary(logBesselI0Scaled(x = x, splineApprox = TRUE) - logBesselI0Scaled(x = x, splineApprox = FALSE)) # Accuracy a1Inv y <- seq(0, 1 - 1e-4, l = 1e3) summary(a1Inv(x = y, splineApprox = TRUE) - a1Inv(x = y, splineApprox = FALSE))# Data employed for log besselI0 scaled x1 <- c(seq(0, 1, by = 1e-4), seq(1 + 1e-2, 10, by = 1e-3), seq(10 + 1e-1, 100, by = 1e-2), seq(100 + 1e0, 1e3, by = 1e0), seq(1000 + 1e1, 5e4, by = 2e1)) logBesselI0ScaledEvalGrid <- log(besselI(x = x1, nu = 0, expon.scaled = TRUE)) # save(list = "logBesselI0ScaledEvalGrid", # file = "logBesselI0ScaledEvalGrid.rda", compress = TRUE) # Data employed for A1 inverse x2 <- rev(c(seq(1e-04, 0.9 - 1e-4, by = 1e-4), seq(0.9, 1 - 1e-05, by = 1e-5))) a1InvEvalGrid <- sapply(x2, function(k) { uniroot(f = function(x) k - besselI(x, nu = 1, expon.scaled = TRUE) / besselI(x, nu = 0, expon.scaled = TRUE), lower = 1e-06, upper = 1e+05, tol = 1e-15)$root }) # save(list = "a1InvEvalGrid", file = "a1InvEvalGrid.rda", compress = TRUE) # Accuracy logBesselI0Scaled x <- seq(0, 1e3, l = 1e3) summary(logBesselI0Scaled(x = x, splineApprox = TRUE) - logBesselI0Scaled(x = x, splineApprox = FALSE)) # Accuracy a1Inv y <- seq(0, 1 - 1e-4, l = 1e3) summary(a1Inv(x = y, splineApprox = TRUE) - a1Inv(x = y, splineApprox = FALSE))
Computation of the loglikelihood for a WN diffusion (with diagonal diffusion matrix) from a sample of initial and final pairs of angles.
logLikWouPairs(x, t, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)logLikWouPairs(x, t, alpha, mu, sigma, rho = 0, maxK = 2L, expTrc = 30)
x |
a matrix of dimension |
t |
either a scalar or a vector of length |
alpha |
vector of length |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
A negative penalty is added if positive definiteness is violated. If the output value is Inf, -100 * N is returned instead.
A scalar giving the final loglikelihood, defined as the sum of the loglikelihood of the initial angles according to the stationary density and the loglikelihood of the transitions from initial to final angles.
set.seed(345567) x <- toPiInt(matrix(rnorm(200, mean = pi), ncol = 4, nrow = 50)) alpha <- c(2, 1, -0.5) mu <- c(0, pi) sigma <- sqrt(c(2, 1)) # The same logLikWouPairs(x = x, t = 0.5, alpha = alpha, mu = mu, sigma = sigma) sum( log(dStatWn2D(x = x[, 1:2], alpha = alpha, mu = mu, sigma = sigma)) + log(dTpdWou2D(x = x[, 3:4], x0 = x[, 1:2], t = 0.5, alpha = alpha, mu = mu, sigma = sigma)) ) # Different times logLikWouPairs(x = x, t = (1:50) / 50, alpha = alpha, mu = mu, sigma = sigma)set.seed(345567) x <- toPiInt(matrix(rnorm(200, mean = pi), ncol = 4, nrow = 50)) alpha <- c(2, 1, -0.5) mu <- c(0, pi) sigma <- sqrt(c(2, 1)) # The same logLikWouPairs(x = x, t = 0.5, alpha = alpha, mu = mu, sigma = sigma) sum( log(dStatWn2D(x = x[, 1:2], alpha = alpha, mu = mu, sigma = sigma)) + log(dTpdWou2D(x = x[, 3:4], x0 = x[, 1:2], t = 0.5, alpha = alpha, mu = mu, sigma = sigma)) ) # Different times logLikWouPairs(x = x, t = (1:50) / 50, alpha = alpha, mu = mu, sigma = sigma)
Computation of the maximum likelihood estimator of the
parameters of the multivariate Ornstein–Uhlenbeck (OU) diffusion
from a discretized trajectory
. The objective
function to minimize is
mleMou(data, delta, alpha = rep(NA, 3), mu = rep(NA, 2), sigma = rep(NA, 2), start, lower = c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01), upper = c(25, 25, 25, pi, pi, 25, 25), ...)mleMou(data, delta, alpha = rep(NA, 3), mu = rep(NA, 2), sigma = rep(NA, 2), start, lower = c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01), upper = c(25, 25, 25, pi, pi, 25, 25), ...)
data |
a matrix of size |
delta |
time discretization step. |
alpha, mu, sigma
|
arguments to fix a parameter to a given value and
perform the estimation on the rest. Defaults to |
start |
starting values, a matrix with |
lower, upper
|
bound for box constraints as in method |
... |
further arguments to be passed to |
The first row in data is not taken into account for
estimation. See mleOu for the univariate case (more efficient).
mleMou only handles p = 2 currently. It imposes that
Sigma is diagonal and handles the parametrization of A by
alphaToA.
Output from mleOptimWrapper.
set.seed(345678) data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = c(1, 1, 0.5), sigma = 1:2), mu = c(1, 1), Sigma = diag((1:2)^2), N = 200, delta = 0.5) mleMou(data = data, delta = 0.5, start = c(1, 1, 0, 1, 1, 1, 2), lower = c(0.1, 0.1, -25, -10, -10, 0.1, 0.1), upper = c(25, 25, 25, 10, 10, 25, 25), maxit = 500) # Fixed sigma and mu mleMou(data = data, delta = 0.5, mu = c(1, 1), sigma = 1:2, start = c(1, 1, 0), lower = c(0.1, 0.1, -25), upper = c(25, 25, 25))set.seed(345678) data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = c(1, 1, 0.5), sigma = 1:2), mu = c(1, 1), Sigma = diag((1:2)^2), N = 200, delta = 0.5) mleMou(data = data, delta = 0.5, start = c(1, 1, 0, 1, 1, 1, 2), lower = c(0.1, 0.1, -25, -10, -10, 0.1, 0.1), upper = c(25, 25, 25, 10, 10, 25, 25), maxit = 500) # Fixed sigma and mu mleMou(data = data, delta = 0.5, mu = c(1, 1), sigma = 1:2, start = c(1, 1, 0), lower = c(0.1, 0.1, -25), upper = c(25, 25, 25))
A convenient wrapper to perform local optimization of the
likelihood function via nlm and optim including several
practical utilities.
mleOptimWrapper(minusLogLik, region = function(pars) list(pars = pars, penalty = 0), penalty = 1e+10, optMethod = "Nelder-Mead", start, lower = rep(-Inf, ncol(start)), upper = rep(Inf, ncol(start)), selectSolution = "lowestLocMin", checkCircular = TRUE, maxit = 500, tol = 1e-05, verbose = 0, eigTol = 1e-04, condTol = 10000, ...)mleOptimWrapper(minusLogLik, region = function(pars) list(pars = pars, penalty = 0), penalty = 1e+10, optMethod = "Nelder-Mead", start, lower = rep(-Inf, ncol(start)), upper = rep(Inf, ncol(start)), selectSolution = "lowestLocMin", checkCircular = TRUE, maxit = 500, tol = 1e-05, verbose = 0, eigTol = 1e-04, condTol = 10000, ...)
minusLogLik |
function computing the minus log-likelihood function.
Must have a single argument containing a vector of length |
region |
function to impose a feasibility region via a penalty. See details. |
penalty |
imposed penalty if value is not finite. |
optMethod |
one of the following strings: |
start |
starting values, a matrix with |
lower, upper
|
bound for box constraints as in method |
selectSolution |
which criterion is used for selecting a solution
among possible ones, either |
checkCircular |
logical indicating whether to automatically treat the
variables with |
maxit |
maximum number of iterations. |
tol |
tolerance for convergence (passed to |
verbose |
an integer from |
eigTol, condTol
|
eigenvalue and condition number tolerance for the
Hessian in order to guarantee a local minimum. Used only if
|
... |
further arguments passed to the |
If checkCircular = TRUE, then the corresponding lower
and upper entries of the circular parameters are set to -Inf
and Inf, respectively, and minusLogLik is called with the
principal value of the circular argument.
If no solution is found satisfying the criterion in selectSolution,
NAs are returned in the elements of the main solution.
The Hessian is only computed if selectSolution = "lowestLocMin".
Region feasibility can be imposed by a function with the same arguments as
minusLogLik that resets pars in to the boundary of the
feasibility region and adds a penalty proportional to the violation of the
feasibility region. Note that this is not the best procedure at all
to solve the constrained optimization problem, but just a relatively
flexible and quick approach (for a more advanced treatment of restrictions,
see
optimization-focused packages). The value must be a list with objects
pars and penalty. By default no region is imposed, i.e.,
region = function(pars) list("pars" = pars, "penalty" = 0). Note that
the Hessian is computed from the unconstrained problem, hence
localMinimumGuaranteed might be FALSE even if a local minimum
to the constrained problem was found.
A list containing the following elements:
par: estimated minimizing parameters
value: value of minusLogLik at the minimum.
convergence: if the optimizer has converged or not.
message: a character string giving any additional information
returned by the optimizer.
eigHessian: eigenvalues of the Hessian at the minimum. Recall
that for the same solution slightly different outputs may be obtained
according to the different computations of the Hessian of nlm and
optim.
localMinimumGuaranteed: tests if the Hessian is positive
definite (all eigenvalues larger than the tolerance eigTol and
condition number smaller than condTol).
solutionsOutput: a list containing the complete output of
the selected method for the different starting values. It includes the
extra objects convergence and localMinimumGuaranteed.
# No local minimum head(mleOptimWrapper(minusLogLik = function(x) -sum((x - 1:4)^2), start = rbind(10:13, 1:2), selectSolution = "lowest")) head(mleOptimWrapper(minusLogLik = function(x) -sum((x - 1:4)^2), start = rbind(10:13, 1:2), selectSolution = "lowestConv")) head(mleOptimWrapper(minusLogLik = function(x) -sum((x - 1:4)^2), start = rbind(10:13, 1:2), selectSolution = "lowestLocMin")) # Local minimum head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:4)^2), start = rbind(10:13), optMethod = "BFGS")) head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:4)^2), start = rbind(10:13), optMethod = "Nelder-Mead")) # Function with several local minimum and a 'spurious' one f <- function(x) 0.75 * (x[1] - 1)^2 - 10 / (0.1 + 0.1 * ((x[1] - 15)^2 + (x[2] - 2)^2)) - 9.5 / (0.1 + 0.1 * ((x[1] - 15)^2 + (x[2] + 2)^2)) plotSurface2D(x = seq(0, 20, l = 100), y = seq(-3, 3, l = 100), f = f) head(mleOptimWrapper(minusLogLik = f, start = rbind(c(15, 2), c(15, -2), c(5, 0)), selectSolution = "lowest")) head(mleOptimWrapper(minusLogLik = f, start = rbind(c(15, 2), c(15, -2), c(5, 0)), selectSolution = "lowestConv")) head(mleOptimWrapper(minusLogLik = f, start = rbind(c(15, 2), c(15, -2), c(5, 0)), selectSolution = "lowestLocMin", eigTol = 0.01)) # With constraint region head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:2)^2), start = rbind(10:11), region = function(pars) { x <- pars[1] y <- pars[2] if (y <= x^2) { return(list("pars" = pars, "penalty" = 0)) } else { return(list("pars" = c(sqrt(y), y), "penalty" = y - x^2)) } }, lower = c(0.5, 1), upper = c(Inf, Inf), optMethod = "Nelder-Mead", selectSolution = "lowest")) head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:2)^2), start = rbind(10:11), lower = c(0.5, 1), upper = c(Inf, Inf),optMethod = "Nelder-Mead"))# No local minimum head(mleOptimWrapper(minusLogLik = function(x) -sum((x - 1:4)^2), start = rbind(10:13, 1:2), selectSolution = "lowest")) head(mleOptimWrapper(minusLogLik = function(x) -sum((x - 1:4)^2), start = rbind(10:13, 1:2), selectSolution = "lowestConv")) head(mleOptimWrapper(minusLogLik = function(x) -sum((x - 1:4)^2), start = rbind(10:13, 1:2), selectSolution = "lowestLocMin")) # Local minimum head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:4)^2), start = rbind(10:13), optMethod = "BFGS")) head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:4)^2), start = rbind(10:13), optMethod = "Nelder-Mead")) # Function with several local minimum and a 'spurious' one f <- function(x) 0.75 * (x[1] - 1)^2 - 10 / (0.1 + 0.1 * ((x[1] - 15)^2 + (x[2] - 2)^2)) - 9.5 / (0.1 + 0.1 * ((x[1] - 15)^2 + (x[2] + 2)^2)) plotSurface2D(x = seq(0, 20, l = 100), y = seq(-3, 3, l = 100), f = f) head(mleOptimWrapper(minusLogLik = f, start = rbind(c(15, 2), c(15, -2), c(5, 0)), selectSolution = "lowest")) head(mleOptimWrapper(minusLogLik = f, start = rbind(c(15, 2), c(15, -2), c(5, 0)), selectSolution = "lowestConv")) head(mleOptimWrapper(minusLogLik = f, start = rbind(c(15, 2), c(15, -2), c(5, 0)), selectSolution = "lowestLocMin", eigTol = 0.01)) # With constraint region head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:2)^2), start = rbind(10:11), region = function(pars) { x <- pars[1] y <- pars[2] if (y <= x^2) { return(list("pars" = pars, "penalty" = 0)) } else { return(list("pars" = c(sqrt(y), y), "penalty" = y - x^2)) } }, lower = c(0.5, 1), upper = c(Inf, Inf), optMethod = "Nelder-Mead", selectSolution = "lowest")) head(mleOptimWrapper(minusLogLik = function(x) sum((x - 1:2)^2), start = rbind(10:11), lower = c(0.5, 1), upper = c(Inf, Inf),optMethod = "Nelder-Mead"))
Computation of the maximum likelihood estimator of the
parameters of the univariate Ornstein–Uhlenbeck (OU) diffusion
from a discretized trajectory
. The objective
function to minimize is
mleOu(data, delta, alpha = NA, mu = NA, sigma = NA, start, lower = c(0.01, -5, 0.01), upper = c(25, 5, 25), ...)mleOu(data, delta, alpha = NA, mu = NA, sigma = NA, start, lower = c(0.01, -5, 0.01), upper = c(25, 5, 25), ...)
data |
a vector of size |
delta |
time discretization step. |
alpha, mu, sigma
|
arguments to fix a parameter to a given value and
perform the estimation on the rest. Defaults to |
start |
starting values, a matrix with |
lower, upper
|
bound for box constraints as in method |
... |
further arguments to be passed to |
The first element in data is not taken into account for
estimation. See mleMou for the multivariate case (less
efficient for dimension one).
Output from mleOptimWrapper.
set.seed(345678) data <- rTrajOu(x0 = 0, alpha = 1, mu = 0, sigma = 1, N = 100, delta = 0.1) mleOu(data = data, delta = 0.1, start = c(2, 1, 2), lower = c(0.1, -10, 0.1), upper = c(25, 10, 25)) # Fixed sigma and mu mleOu(data = data, delta = 0.1, mu = 0, sigma = 1, start = 2, lower = 0.1, upper = 25, optMethod = "nlm")set.seed(345678) data <- rTrajOu(x0 = 0, alpha = 1, mu = 0, sigma = 1, N = 100, delta = 0.1) mleOu(data = data, delta = 0.1, start = c(2, 1, 2), lower = c(0.1, -10, 0.1), upper = c(25, 10, 25)) # Fixed sigma and mu mleOu(data = data, delta = 0.1, mu = 0, sigma = 1, start = 2, lower = 0.1, upper = 25, optMethod = "nlm")
Maximum Likelihood Estimation (MLE) for arbitrary diffusions, based on the transition probability density (tpd) obtained as the numerical solution of the Fokker–Planck Partial Differential Equation (PDE) in 1D.
mlePde1D(data, delta, b, sigma2, Mx = 500, Mt = ceiling(100 * delta), sdInitial = 0.1, linearBinning = FALSE, start, lower, upper, ...)mlePde1D(data, delta, b, sigma2, Mx = 500, Mt = ceiling(100 * delta), sdInitial = 0.1, linearBinning = FALSE, start, lower, upper, ...)
data |
a vector of size |
delta |
time discretization step. |
b |
drift function. Must return a vector of the same size as its argument. |
sigma2 |
function giving the squared diffusion coefficient. Must return a vector of the same size as its argument. |
Mx |
size of the equispaced spatial grid in |
Mt |
size of the time grid in |
sdInitial |
the standard deviation of the concentrated WN giving the initial condition. |
linearBinning |
flag to indicate whether linear binning should be applied for the initial values of the tpd, instead of usual simple binning (cheaper). Linear binning is always done in the evaluation of the tpd. |
start |
starting values, a matrix with |
lower, upper
|
bound for box constraints as in method |
... |
Further parameters passed to |
See Sections 3.4.1 and 3.4.4 in García-Portugués et al. (2019) for details.
Output from mleOptimWrapper.
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
# Test in OU alpha <- 2 mu <- 0 sigma <- 1 set.seed(234567) traj <- rTrajOu(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 500, delta = 0.5) b <- function(x, pars) pars[1] * (pars[2] - x) sigma2 <- function(x, pars) rep(pars[3]^2, length(x)) exactOu <- mleOu(traj, delta = 0.5, start = c(1, 1, 2), lower = c(0.1, -pi, 0.1), upper = c(10, pi, 10)) pdeOu <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = b, sigma2 = sigma2, start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), verbose = 2) pdeOuLin <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = b, sigma2 = sigma2, start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), linearBinning = TRUE, verbose = 2) head(exactOu) head(pdeOu) head(pdeOuLin) # Test in WN diffusion alpha <- 2 mu <- 0 sigma <- 1 set.seed(234567) traj <- rTrajWn1D(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 500, delta = 0.5) exactOu <- mleOu(traj, delta = 0.5, start = c(1, 1, 2), lower = c(0.1, -pi, 0.1), upper = c(10, pi, 10)) pdeWn <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2], sigma = pars[3]), sigma2 = function(x, pars) rep(pars[3]^2, length(x)), start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), verbose = 2) pdeWnLin <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2], sigma = pars[3]), sigma2 = function(x, pars) rep(pars[3]^2, length(x)), start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), linearBinning = TRUE, verbose = 2) head(exactOu) head(pdeWn) head(pdeWnLin)# Test in OU alpha <- 2 mu <- 0 sigma <- 1 set.seed(234567) traj <- rTrajOu(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 500, delta = 0.5) b <- function(x, pars) pars[1] * (pars[2] - x) sigma2 <- function(x, pars) rep(pars[3]^2, length(x)) exactOu <- mleOu(traj, delta = 0.5, start = c(1, 1, 2), lower = c(0.1, -pi, 0.1), upper = c(10, pi, 10)) pdeOu <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = b, sigma2 = sigma2, start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), verbose = 2) pdeOuLin <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = b, sigma2 = sigma2, start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), linearBinning = TRUE, verbose = 2) head(exactOu) head(pdeOu) head(pdeOuLin) # Test in WN diffusion alpha <- 2 mu <- 0 sigma <- 1 set.seed(234567) traj <- rTrajWn1D(x0 = 0, alpha = alpha, mu = mu, sigma = sigma, N = 500, delta = 0.5) exactOu <- mleOu(traj, delta = 0.5, start = c(1, 1, 2), lower = c(0.1, -pi, 0.1), upper = c(10, pi, 10)) pdeWn <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2], sigma = pars[3]), sigma2 = function(x, pars) rep(pars[3]^2, length(x)), start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), verbose = 2) pdeWnLin <- mlePde1D(data = traj, delta = 0.5, Mx = 100, Mt = 100, b = function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2], sigma = pars[3]), sigma2 = function(x, pars) rep(pars[3]^2, length(x)), start = c(1, 1, 2), lower = c(0.1, -pi, -10), upper = c(10, pi, 10), linearBinning = TRUE, verbose = 2) head(exactOu) head(pdeWn) head(pdeWnLin)
Maximum Likelihood Estimation (MLE) for arbitrary diffusions, based on the transition probability density (tpd) obtained as the numerical solution of the Fokker–Planck Partial Differential Equation (PDE) in 2D.
mlePde2D(data, delta, b, sigma2, Mx = 50, My = 50, Mt = ceiling(100 * delta), sdInitial = 0.1, linearBinning = FALSE, start, lower, upper, ...)mlePde2D(data, delta, b, sigma2, Mx = 50, My = 50, Mt = ceiling(100 * delta), sdInitial = 0.1, linearBinning = FALSE, start, lower, upper, ...)
data |
a matrix of dimension |
delta |
discretization step. |
b |
drift function. Must return a vector of the same size as its argument. |
sigma2 |
function giving the diagonal of the diffusion matrix. Must return a vector of the same size as its argument. |
Mx, My
|
sizes of the equispaced spatial grids in |
Mt |
size of the time grid in |
sdInitial |
standard deviations of the concentrated WN giving the initial condition. |
linearBinning |
flag to indicate whether linear binning should be applied for the initial values of the tpd, instead of usual simple binning (cheaper). Linear binning is always done in the evaluation of the tpd. |
start |
starting values, a matrix with |
lower, upper
|
bound for box constraints as in method |
... |
further parameters passed to |
See Sections 3.4.2 and 3.4.4 in García-Portugués et al. (2019) for
details. The function currently includes the region function for
imposing a feasibility region on the parameters of the bivariate WN
diffusion.
Output from mleOptimWrapper.
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
# Test in OU process alpha <- c(1, 2, -0.5) mu <- c(0, 0) sigma <- c(0.5, 1) set.seed(2334567) data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = alpha, sigma = sigma), mu = mu, Sigma = diag(sigma^2), N = 500, delta = 0.5) b <- function(x, pars) sweep(-x, 2, pars[4:5], "+") %*% t(alphaToA(alpha = pars[1:3], sigma = sigma)) sigma2 <- function(x, pars) repRow(sigma^2, nrow(x)) exactOu <- mleMou(data = data, delta = 0.5, sigma = sigma, start = c(1, 1, 0, 2, 2), lower = c(0.1, 0.1, -25, -10, -10), upper = c(25, 25, 25, 10, 10)) head(exactOu, 2) pdeOu <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 10, My = 10, Mt = 10, start = rbind(c(1, 1, 0, 2, 2)), lower = c(0.1, 0.1, -25, -10, -10), upper = c(25, 25, 25, 10, 10), verbose = 2) head(pdeOu, 2) pdeOuLin <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 10, My = 10, Mt = 10, start = rbind(c(1, 1, 0, 2, 2)), lower = c(0.1, 0.1, -25, -10, -10), upper = c(25, 25, 25, 10, 10), verbose = 2, linearBinning = TRUE) head(pdeOuLin, 2) # Test in WN diffusion alpha <- c(1, 0.5, 0.25) mu <- c(0, 0) sigma <- c(2, 1) set.seed(234567) data <- rTrajWn2D(x0 = c(0, 0), alpha = alpha, mu = mu, sigma = sigma, N = 200, delta = 0.5) b <- function(x, pars) driftWn2D(x = x, A = alphaToA(alpha = pars[1:3], sigma = sigma), mu = pars[4:5], sigma = sigma) sigma2 <- function(x, pars) repRow(sigma^2, nrow(x)) exactOu <- mleMou(data = data, delta = 0.5, sigma = sigma, start = c(1, 1, 0, 1, 1), lower = c(0.1, 0.1, -25, -25, -25), upper = c(25, 25, 25, 25, 25), optMethod = "nlm") pdeWn <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 20, My = 20, Mt = 10, start = rbind(c(1, 1, 0, 1, 1)), lower = c(0.1, 0.1, -25, -25, -25), upper = c(25, 25, 25, 25, 25), verbose = 2, optMethod = "nlm") pdeWnLin <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 20, My = 20, Mt = 10, start = rbind(c(1, 1, 0, 1, 1)), lower = c(0.1, 0.1, -25, -25, -25), upper = c(25, 25, 25, 25, 25), verbose = 2, linearBinning = TRUE) head(exactOu) head(pdeOu) head(pdeOuLin)# Test in OU process alpha <- c(1, 2, -0.5) mu <- c(0, 0) sigma <- c(0.5, 1) set.seed(2334567) data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = alpha, sigma = sigma), mu = mu, Sigma = diag(sigma^2), N = 500, delta = 0.5) b <- function(x, pars) sweep(-x, 2, pars[4:5], "+") %*% t(alphaToA(alpha = pars[1:3], sigma = sigma)) sigma2 <- function(x, pars) repRow(sigma^2, nrow(x)) exactOu <- mleMou(data = data, delta = 0.5, sigma = sigma, start = c(1, 1, 0, 2, 2), lower = c(0.1, 0.1, -25, -10, -10), upper = c(25, 25, 25, 10, 10)) head(exactOu, 2) pdeOu <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 10, My = 10, Mt = 10, start = rbind(c(1, 1, 0, 2, 2)), lower = c(0.1, 0.1, -25, -10, -10), upper = c(25, 25, 25, 10, 10), verbose = 2) head(pdeOu, 2) pdeOuLin <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 10, My = 10, Mt = 10, start = rbind(c(1, 1, 0, 2, 2)), lower = c(0.1, 0.1, -25, -10, -10), upper = c(25, 25, 25, 10, 10), verbose = 2, linearBinning = TRUE) head(pdeOuLin, 2) # Test in WN diffusion alpha <- c(1, 0.5, 0.25) mu <- c(0, 0) sigma <- c(2, 1) set.seed(234567) data <- rTrajWn2D(x0 = c(0, 0), alpha = alpha, mu = mu, sigma = sigma, N = 200, delta = 0.5) b <- function(x, pars) driftWn2D(x = x, A = alphaToA(alpha = pars[1:3], sigma = sigma), mu = pars[4:5], sigma = sigma) sigma2 <- function(x, pars) repRow(sigma^2, nrow(x)) exactOu <- mleMou(data = data, delta = 0.5, sigma = sigma, start = c(1, 1, 0, 1, 1), lower = c(0.1, 0.1, -25, -25, -25), upper = c(25, 25, 25, 25, 25), optMethod = "nlm") pdeWn <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 20, My = 20, Mt = 10, start = rbind(c(1, 1, 0, 1, 1)), lower = c(0.1, 0.1, -25, -25, -25), upper = c(25, 25, 25, 25, 25), verbose = 2, optMethod = "nlm") pdeWnLin <- mlePde2D(data = data, delta = 0.5, b = b, sigma2 = sigma2, Mx = 20, My = 20, Mt = 10, start = rbind(c(1, 1, 0, 1, 1)), lower = c(0.1, 0.1, -25, -25, -25), upper = c(25, 25, 25, 25, 25), verbose = 2, linearBinning = TRUE) head(exactOu) head(pdeOu) head(pdeOuLin)
Quadrature rules for definite integrals over intervals in 1D,
, rectangles in 2D, and cubes in 3D,
.
The trapezoidal rules assume that the function is periodic, whereas the
Simpson rules work for arbitrary functions.
periodicTrapRule1D(fx, endsMatch = FALSE, na.rm = TRUE, lengthInterval = 2 * pi) periodicTrapRule2D(fxy, endsMatch = FALSE, na.rm = TRUE, lengthInterval = rep(2 * pi, 2)) periodicTrapRule3D(fxyz, endsMatch = FALSE, na.rm = TRUE, lengthInterval = rep(2 * pi, 3)) integrateSimp1D(fx, lengthInterval = 2 * pi, na.rm = TRUE) integrateSimp2D(fxy, lengthInterval = rep(2 * pi, 2), na.rm = TRUE) integrateSimp3D(fxyz, lengthInterval = rep(2 * pi, 3), na.rm = TRUE)periodicTrapRule1D(fx, endsMatch = FALSE, na.rm = TRUE, lengthInterval = 2 * pi) periodicTrapRule2D(fxy, endsMatch = FALSE, na.rm = TRUE, lengthInterval = rep(2 * pi, 2)) periodicTrapRule3D(fxyz, endsMatch = FALSE, na.rm = TRUE, lengthInterval = rep(2 * pi, 3)) integrateSimp1D(fx, lengthInterval = 2 * pi, na.rm = TRUE) integrateSimp2D(fxy, lengthInterval = rep(2 * pi, 2), na.rm = TRUE) integrateSimp3D(fxyz, lengthInterval = rep(2 * pi, 3), na.rm = TRUE)
fx |
vector containing the evaluation of the function to integrate over
a uniform grid in |
endsMatch |
flag to indicate whether the values of the last entries of
|
na.rm |
logical. Should missing values (including |
lengthInterval |
vector containing the lengths of the intervals of integration. |
fxy |
matrix containing the evaluation of the function to integrate
over a uniform grid in |
fxyz |
three dimensional array containing the evaluation of the
function to integrate over a uniform grid in
|
The simple trapezoidal rule has a very good performance for periodic functions in 1D and 2D(order of error ). The higher dimensional extensions are obtained by iterative usage of the 1D rules.
The value of the integral.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. (1996). Numerical Recipes in Fortran 77: The Art of Scientific Computing (Vol. 1 of Fortran Numerical Recipes). Cambridge University Press, Cambridge.
# In 1D. True value: 3.55099937 N <- 21 grid <- seq(-pi, pi, l = N) fx <- sin(grid)^2 * exp(cos(grid)) periodicTrapRule1D(fx = fx, endsMatch = TRUE) periodicTrapRule1D(fx = fx[-N], endsMatch = FALSE) integrateSimp1D(fx = fx, lengthInterval = 2 * pi) integrateSimp1D(fx = fx[-N]) # Worse, of course # In 2D. True value: 22.31159 fxy <- outer(grid, grid, function(x, y) (sin(x)^2 * exp(cos(x)) + sin(y)^2 * exp(cos(y))) / 2) periodicTrapRule2D(fxy = fxy, endsMatch = TRUE) periodicTrapRule2D(fxy = fxy[-N, -N], endsMatch = FALSE) periodicTrapRule1D(apply(fxy[-N, -N], 1, periodicTrapRule1D)) integrateSimp2D(fxy = fxy) integrateSimp1D(apply(fxy, 1, integrateSimp1D)) # In 3D. True value: 140.1878 fxyz <- array(fxy, dim = c(N, N, N)) for (i in 1:N) fxyz[i, , ] <- fxy periodicTrapRule3D(fxyz = fxyz, endsMatch = TRUE) integrateSimp3D(fxyz = fxyz)# In 1D. True value: 3.55099937 N <- 21 grid <- seq(-pi, pi, l = N) fx <- sin(grid)^2 * exp(cos(grid)) periodicTrapRule1D(fx = fx, endsMatch = TRUE) periodicTrapRule1D(fx = fx[-N], endsMatch = FALSE) integrateSimp1D(fx = fx, lengthInterval = 2 * pi) integrateSimp1D(fx = fx[-N]) # Worse, of course # In 2D. True value: 22.31159 fxy <- outer(grid, grid, function(x, y) (sin(x)^2 * exp(cos(x)) + sin(y)^2 * exp(cos(y))) / 2) periodicTrapRule2D(fxy = fxy, endsMatch = TRUE) periodicTrapRule2D(fxy = fxy[-N, -N], endsMatch = FALSE) periodicTrapRule1D(apply(fxy[-N, -N], 1, periodicTrapRule1D)) integrateSimp2D(fxy = fxy) integrateSimp1D(apply(fxy, 1, integrateSimp1D)) # In 3D. True value: 140.1878 fxyz <- array(fxy, dim = c(N, N, N)) for (i in 1:N) fxyz[i, , ] <- fxy periodicTrapRule3D(fxyz = fxyz, endsMatch = TRUE) integrateSimp3D(fxyz = fxyz)
Maximum pseudo-likelihood using the Euler and Shoji–Ozaki pseudo-likelihoods.
psMle(data, delta, method = c("E", "SO", "SO2"), b, jac.b, sigma2, b1, b2, start, lower, upper, circular = TRUE, maxK = 2, vmApprox = FALSE, ...)psMle(data, delta, method = c("E", "SO", "SO2"), b, jac.b, sigma2, b1, b2, start, lower, upper, circular = TRUE, maxK = 2, vmApprox = FALSE, ...)
data |
a matrix of dimension |
delta |
discretization step. |
method |
a string for choosing |
b |
drift function. Must return a matrix of the same size as |
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 |
b1 |
first derivative of the drift function (univariate). Must return
a vector of the same length as |
b2 |
second derivative of the drift function (univariate). Must return
a vector of the same length as |
start |
starting values, a matrix with |
lower, upper
|
bound for box constraints as in method |
circular |
flag to indicate circular data. |
maxK |
maximum absolute winding number used if |
vmApprox |
flag to indicate von Mises approximation to wrapped normal.
See |
... |
further parameters passed to |
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.
Output from mleOptimWrapper.
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
# Example in 1D delta <- 0.5 pars <- c(0.25, 0, 2) set.seed(12345678) samp <- rTrajWn1D(x0 = 0, alpha = pars[1], mu = pars[2], sigma = pars[3], N = 100, delta = delta) b <- function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2], sigma = pars[3], maxK = 2, expTrc = 30) b1 <- function(x, pars, h = 1e-4) { l <- length(x) res <- b(x = c(x + h, x - h), pars = pars) drop(res[1:l] - res[(l + 1):(2 * l)])/(2 * h) } b2 <- function(x, pars, h = 1e-4) { l <- length(x) res <- b(x = c(x + h, x, x - h), pars = pars) drop(res[1:l] - 2 * res[(l + 1):(2 * l)] + res[(2 * l + 1):(3 * l)])/(h^2) } sigma2 <- function(x, pars) rep(pars[3]^2, length(x)) lower <- c(0.1, -pi, 0.1) upper <- c(10, pi, 10) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper, vmApprox = TRUE) psMle(data = samp, delta = delta, method = "SO2", b = b, b1 = b1, b2 = b2, sigma2 = sigma2, start = pars, lower = lower, upper = upper) psMle(data = samp, delta = delta, method = "SO2", b = b, b1 = b1, b2 = b2, sigma2 = sigma2, start = pars, lower = lower, upper = upper, vmApprox = TRUE) psMle(data = samp, delta = delta, method = "SO", b = b, b1 = b1, lower = lower, upper = upper, sigma2 = sigma2, start = pars) approxMleWn1D(data = samp, delta = delta, start = pars) mlePde1D(data = samp, delta = delta, b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) # Example in 2D delta <- 0.5 pars <- c(1, 0.5, 0, 0, 0, 1, 2) set.seed(12345678) samp <- rTrajWn2D(x0 = c(0, 0), alpha = pars[1:3], mu = pars[4:5], sigma = pars[6:7], N = 100, delta = delta) b <- function(x, pars) driftWn2D(x = x, A = alphaToA(alpha = pars[1:3], sigma = pars[6:7]), mu = pars[4:5], sigma = pars[6:7], maxK = 2, expTrc = 30) jac.b <- function(x, pars, h = 1e-4) { l <- nrow(x) res <- b(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)), pars = pars) cbind(res[1:l, ] - res[(l + 1):(2 * l), ], res[2 * l + 1:l, ] - res[2 * l + (l + 1):(2 * l), ]) / (2 * h) } sigma2 <- function(x, pars) matrix(pars[6:7]^2, nrow = length(x) / 2L, ncol = 2) lower <- c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01) upper <- c(25, 25, 25, pi, pi, 25, 25) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper, vmApprox = TRUE) psMle(data = samp, delta = delta, method = "SO", b = b, jac.b = jac.b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) approxMleWn2D(data = samp, delta = delta, start = c(pars, 0)) # Set maxit = 5 to test and avoid a very long evaluation mlePde2D(data = samp, delta = delta, b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper, maxit = 5)# Example in 1D delta <- 0.5 pars <- c(0.25, 0, 2) set.seed(12345678) samp <- rTrajWn1D(x0 = 0, alpha = pars[1], mu = pars[2], sigma = pars[3], N = 100, delta = delta) b <- function(x, pars) driftWn1D(x = x, alpha = pars[1], mu = pars[2], sigma = pars[3], maxK = 2, expTrc = 30) b1 <- function(x, pars, h = 1e-4) { l <- length(x) res <- b(x = c(x + h, x - h), pars = pars) drop(res[1:l] - res[(l + 1):(2 * l)])/(2 * h) } b2 <- function(x, pars, h = 1e-4) { l <- length(x) res <- b(x = c(x + h, x, x - h), pars = pars) drop(res[1:l] - 2 * res[(l + 1):(2 * l)] + res[(2 * l + 1):(3 * l)])/(h^2) } sigma2 <- function(x, pars) rep(pars[3]^2, length(x)) lower <- c(0.1, -pi, 0.1) upper <- c(10, pi, 10) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper, vmApprox = TRUE) psMle(data = samp, delta = delta, method = "SO2", b = b, b1 = b1, b2 = b2, sigma2 = sigma2, start = pars, lower = lower, upper = upper) psMle(data = samp, delta = delta, method = "SO2", b = b, b1 = b1, b2 = b2, sigma2 = sigma2, start = pars, lower = lower, upper = upper, vmApprox = TRUE) psMle(data = samp, delta = delta, method = "SO", b = b, b1 = b1, lower = lower, upper = upper, sigma2 = sigma2, start = pars) approxMleWn1D(data = samp, delta = delta, start = pars) mlePde1D(data = samp, delta = delta, b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) # Example in 2D delta <- 0.5 pars <- c(1, 0.5, 0, 0, 0, 1, 2) set.seed(12345678) samp <- rTrajWn2D(x0 = c(0, 0), alpha = pars[1:3], mu = pars[4:5], sigma = pars[6:7], N = 100, delta = delta) b <- function(x, pars) driftWn2D(x = x, A = alphaToA(alpha = pars[1:3], sigma = pars[6:7]), mu = pars[4:5], sigma = pars[6:7], maxK = 2, expTrc = 30) jac.b <- function(x, pars, h = 1e-4) { l <- nrow(x) res <- b(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)), pars = pars) cbind(res[1:l, ] - res[(l + 1):(2 * l), ], res[2 * l + 1:l, ] - res[2 * l + (l + 1):(2 * l), ]) / (2 * h) } sigma2 <- function(x, pars) matrix(pars[6:7]^2, nrow = length(x) / 2L, ncol = 2) lower <- c(0.01, 0.01, -25, -pi, -pi, 0.01, 0.01) upper <- c(25, 25, 25, pi, pi, 25, 25) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) psMle(data = samp, delta = delta, method = "E", b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper, vmApprox = TRUE) psMle(data = samp, delta = delta, method = "SO", b = b, jac.b = jac.b, sigma2 = sigma2, start = pars, lower = lower, upper = upper) approxMleWn2D(data = samp, delta = delta, start = c(pars, 0)) # Set maxit = 5 to test and avoid a very long evaluation mlePde2D(data = samp, delta = delta, b = b, sigma2 = sigma2, start = pars, lower = lower, upper = upper, maxit = 5)
Simulates from the stationary density of the WN diffusion in 2D.
rStatWn2D(n, mu, alpha, sigma, rho = 0)rStatWn2D(n, mu, alpha, sigma, rho = 0)
n |
sample size. |
mu |
a vector of length |
alpha |
vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
A matrix of dimension c(n, 2) containing the samples from the stationary distribution.
set.seed(345567) alpha <- c(2, 1, -1) sigma <- c(1.5, 2) Sigma <- diag(sigma^2) A <- alphaToA(alpha = alpha, sigma = sigma) mu <- c(pi, pi) plot(rStatWn2D(n = 1000, mu = mu, alpha = alpha, sigma = sigma)) points(toPiInt(mvtnorm::rmvnorm(n = 1000, mean = mu, sigma = solve(A) %*% Sigma / 2, method = "chol")), col = 2)set.seed(345567) alpha <- c(2, 1, -1) sigma <- c(1.5, 2) Sigma <- diag(sigma^2) A <- alphaToA(alpha = alpha, sigma = sigma) mu <- c(pi, pi) plot(rStatWn2D(n = 1000, mu = mu, alpha = alpha, sigma = sigma)) points(toPiInt(mvtnorm::rmvnorm(n = 1000, mean = mu, sigma = solve(A) %*% Sigma / 2, method = "chol")), col = 2)
Simulates from the approximate transition density of the WN diffusion in 2D.
rTpdWn2D(n, x0, t, mu, alpha, sigma, rho = 0, maxK = 2L, expTrc = 30)rTpdWn2D(n, x0, t, mu, alpha, sigma, rho = 0, maxK = 2L, expTrc = 30)
n |
sample size. |
x0 |
a matrix of dimension |
t |
vector of length |
mu |
a vector of length |
alpha |
vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
An array of dimension c(n, 2, nx0) containing the n samples of the transition distribution,
conditioned on that the process was at x0 at t instants ago. The samples are all in .
alpha <- c(3, 2, -1) sigma <- c(0.5, 1) mu <- c(pi, pi) x <- seq(-pi, pi, l = 100) t <- 0.5 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, mu = mu, alpha = alpha, sigma = sigma, maxK = 2, expTrc = 30), nrow = 100, ncol = 100), zlim = c(0, 0.5)) points(rTpdWn2D(n = 500, x0 = rbind(c(0, 0)), t = t, mu = mu, alpha = alpha, sigma = sigma)[, , 1], col = 3) points(stepAheadWn2D(x0 = rbind(c(0, 0)), delta = t / 500, A = alphaToA(alpha = alpha, sigma = sigma), mu = mu, sigma = sigma, N = 500, M = 500, maxK = 2, expTrc = 30), col = 4)alpha <- c(3, 2, -1) sigma <- c(0.5, 1) mu <- c(pi, pi) x <- seq(-pi, pi, l = 100) t <- 0.5 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, mu = mu, alpha = alpha, sigma = sigma, maxK = 2, expTrc = 30), nrow = 100, ncol = 100), zlim = c(0, 0.5)) points(rTpdWn2D(n = 500, x0 = rbind(c(0, 0)), t = t, mu = mu, alpha = alpha, sigma = sigma)[, , 1], col = 3) points(stepAheadWn2D(x0 = rbind(c(0, 0)), delta = t / 500, A = alphaToA(alpha = alpha, sigma = sigma), mu = mu, sigma = sigma, N = 500, M = 500, maxK = 2, expTrc = 30), col = 4)
Simulation of an arbitrary Langevin diffusion in dimension
p by subsampling a fine trajectory obtained by the Euler
discretization.
rTrajLangevin(x0, drift, SigDif, N = 100, delta = 0.01, NFine = ceiling(N * delta/deltaFine), deltaFine = min(delta/100, 0.001), circular = TRUE, ...)rTrajLangevin(x0, drift, SigDif, N = 100, delta = 0.01, NFine = ceiling(N * delta/deltaFine), deltaFine = min(delta/100, 0.001), circular = TRUE, ...)
x0 |
vector of length |
drift |
drift for the diffusion. |
SigDif |
matrix of size |
N |
number of discretization steps in the resulting trajectory. |
delta |
discretization step. |
NFine |
number of discretization steps for the fine trajectory. Must
be larger than |
deltaFine |
discretization step for the fine trajectory. Must be
smaller than |
circular |
whether to wrap the resulting trajectory to
|
... |
parameters to be passed to |
The fine trajectory is subsampled using the indexes
seq(1, NFine + 1, by = NFine / N).
A vector of length N + 1 containing x0 in the first
entry and the discretized trajectory.
isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { # 1D manipulate::manipulate({ x <- seq(0, N * delta, by = delta) plot(x, x, ylim = c(-pi, pi), type = "n", ylab = expression(X[t]), xlab = "t") linesCirc(x, rTrajLangevin(x0 = 0, drift = driftJp, SigDif = sigma, alpha = alpha, mu = 0, psi = psi, N = N, delta = 0.01)) }, delta = manipulate::slider(0.01, 5.01, step = 0.1), N = manipulate::slider(10, 500, step = 10, initial = 200), alpha = manipulate::slider(0.01, 5, step = 0.1, initial = 1), psi = manipulate::slider(-2, 2, step = 0.1, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.1, initial = 1)) # 2D samp <- rTrajLangevin(x0 = c(0, 0), drift = driftMvm, alpha = c(1, 1), mu = c(2, -1), A = diag(rep(0, 2)), SigDif = diag(rep(1, 2)), N = 1000, delta = 0.1) plot(samp, xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 19, cex = 0.25, xlab = expression(X[t]), ylab = expression(Y[t]), col = rainbow(1000)) linesTorus(samp[, 1], samp[, 2], col = rainbow(1000)) }isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { # 1D manipulate::manipulate({ x <- seq(0, N * delta, by = delta) plot(x, x, ylim = c(-pi, pi), type = "n", ylab = expression(X[t]), xlab = "t") linesCirc(x, rTrajLangevin(x0 = 0, drift = driftJp, SigDif = sigma, alpha = alpha, mu = 0, psi = psi, N = N, delta = 0.01)) }, delta = manipulate::slider(0.01, 5.01, step = 0.1), N = manipulate::slider(10, 500, step = 10, initial = 200), alpha = manipulate::slider(0.01, 5, step = 0.1, initial = 1), psi = manipulate::slider(-2, 2, step = 0.1, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.1, initial = 1)) # 2D samp <- rTrajLangevin(x0 = c(0, 0), drift = driftMvm, alpha = c(1, 1), mu = c(2, -1), A = diag(rep(0, 2)), SigDif = diag(rep(1, 2)), N = 1000, delta = 0.1) plot(samp, xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 19, cex = 0.25, xlab = expression(X[t]), ylab = expression(Y[t]), col = rainbow(1000)) linesTorus(samp[, 1], samp[, 2], col = rainbow(1000)) }
Simulation of trajectories of the multivariate Ornstein–Uhlenbeck (OU) diffusion
using the exact transition probability density.
rTrajMou(x0, A, mu, Sigma, N = 100, delta = 0.001)rTrajMou(x0, A, mu, Sigma, N = 100, delta = 0.001)
x0 |
a vector of length |
A |
the drift matrix, of size |
mu |
unconditional mean of the diffusion, a vector of length |
Sigma |
square of the diffusion matrix, a matrix of size |
N |
number of discretization steps in the resulting trajectory. |
delta |
time discretization step. |
The law of the discretized trajectory at each time step is
a multivariate normal with mean meantMou and covariance
matrix covtMou. See rTrajOu for the univariate
case (more efficient).
solve(A) %*% Sigma has to be a covariance matrix (symmetric and
positive definite) in order to have a proper transition density. For the
bivariate case, this can be ensured with the alphaToA function.
In the multivariate case, it is ensured if Sigma is isotropic and
A is a covariance matrix.
A matrix of size c(N + 1, p) containing x0 in the
first row and the exact discretized trajectory on the remaining rows.
set.seed(987658) data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = c(1, 2, 0.5), sigma = 1:2), mu = c(1, 1), Sigma = diag((1:2)^2), N = 200, delta = 0.1) plot(data, pch = 19, col = rainbow(201), cex = 0.25) arrows(x0 = data[-201, 1], y0 = data[-201, 2], x1 = data[-1, 1], y1 = data[-1, 2], col = rainbow(201), angle = 10, length = 0.1)set.seed(987658) data <- rTrajMou(x0 = c(0, 0), A = alphaToA(alpha = c(1, 2, 0.5), sigma = 1:2), mu = c(1, 1), Sigma = diag((1:2)^2), N = 200, delta = 0.1) plot(data, pch = 19, col = rainbow(201), cex = 0.25) arrows(x0 = data[-201, 1], y0 = data[-201, 2], x1 = data[-1, 1], y1 = data[-1, 2], col = rainbow(201), angle = 10, length = 0.1)
Simulation of trajectories of the univariate Ornstein–Uhlenbeck (OU) diffusion
using the exact transition probability density.
rTrajOu(x0, alpha, mu, sigma, N = 100, delta = 0.001)rTrajOu(x0, alpha, mu, sigma, N = 100, delta = 0.001)
x0 |
initial point. |
alpha |
strength of the drift. |
mu |
unconditional mean of the diffusion. |
sigma |
diffusion coefficient. |
N |
number of discretization steps in the resulting trajectory. |
delta |
time discretization step. |
The law of the discretized trajectory is a multivariate normal
with mean meantOu and covariance matrix covstOu.
See rTrajMou for the multivariate case (less efficient for
dimension one).
A vector of length N + 1 containing x0 in the first
entry and the exact discretized trajectory on the remaining elements.
isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate({ set.seed(345678); plot(seq(0, N * delta, by = delta), rTrajOu(x0 = 0, alpha = alpha, mu = 0, sigma = sigma, N = N, delta = delta), ylim = c(-4, 4), type = "l", ylab = expression(X[t]), xlab = "t") }, delta = manipulate::slider(0.01, 5.01, step = 0.1), N = manipulate::slider(10, 500, step = 10, initial = 200), alpha = manipulate::slider(0.01, 5, step = 0.1, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.1, initial = 1)) }isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate({ set.seed(345678); plot(seq(0, N * delta, by = delta), rTrajOu(x0 = 0, alpha = alpha, mu = 0, sigma = sigma, N = N, delta = delta), ylim = c(-4, 4), type = "l", ylab = expression(X[t]), xlab = "t") }, delta = manipulate::slider(0.01, 5.01, step = 0.1), N = manipulate::slider(10, 500, step = 10, initial = 200), alpha = manipulate::slider(0.01, 5, step = 0.1, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.1, initial = 1)) }
Simulation of the Wrapped Normal (WN) diffusion in 1D by subsampling a fine trajectory obtained by the Euler discretization.
rTrajWn1D(x0, alpha, mu, sigma, N = 100, delta = 0.01, NFine = ceiling(N * delta/deltaFine), deltaFine = min(delta/100, 0.001))rTrajWn1D(x0, alpha, mu, sigma, N = 100, delta = 0.01, NFine = ceiling(N * delta/deltaFine), deltaFine = min(delta/100, 0.001))
x0 |
initial point. |
alpha |
drift parameter. |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
N |
number of discretization steps in the resulting trajectory. |
delta |
discretization step. |
NFine |
number of discretization steps for the fine trajectory. Must
be larger than |
deltaFine |
discretization step for the fine trajectory. Must be
smaller than |
The fine trajectory is subsampled using the indexes
seq(1, NFine + 1, by = NFine / N).
A vector of length N + 1 containing x0 in the first
entry and the discretized trajectory.
isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate({ x <- seq(0, N * delta, by = delta) plot(x, x, ylim = c(-pi, pi), type = "n", ylab = expression(X[t]), xlab = "t") linesCirc(x, rTrajWn1D(x0 = 0, alpha = alpha, mu = 0, sigma = sigma, N = N, delta = 0.01)) }, delta = slider(0.01, 5.01, step = 0.1), N = manipulate::slider(10, 500, step = 10, initial = 200), alpha = manipulate::slider(0.01, 5, step = 0.1, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.1, initial = 1)) }isRStudio <- identical(.Platform$GUI, "RStudio") if (isRStudio) { manipulate::manipulate({ x <- seq(0, N * delta, by = delta) plot(x, x, ylim = c(-pi, pi), type = "n", ylab = expression(X[t]), xlab = "t") linesCirc(x, rTrajWn1D(x0 = 0, alpha = alpha, mu = 0, sigma = sigma, N = N, delta = 0.01)) }, delta = slider(0.01, 5.01, step = 0.1), N = manipulate::slider(10, 500, step = 10, initial = 200), alpha = manipulate::slider(0.01, 5, step = 0.1, initial = 1), sigma = manipulate::slider(0.01, 5, step = 0.1, initial = 1)) }
Simulation of the Wrapped Normal (WN) diffusion in 2D by subsampling a fine trajectory obtained by the Euler discretization.
rTrajWn2D(x0, alpha, mu, sigma, rho = 0, N = 100, delta = 0.01, NFine = ceiling(N * delta/deltaFine), deltaFine = min(delta/100, 0.001))rTrajWn2D(x0, alpha, mu, sigma, rho = 0, N = 100, delta = 0.01, NFine = ceiling(N * delta/deltaFine), deltaFine = min(delta/100, 0.001))
x0 |
vector of length |
alpha |
vector of length |
mu |
a vector of length |
sigma |
vector of length |
rho |
correlation coefficient of |
N |
number of discretization steps in the resulting trajectory. |
delta |
discretization step. |
NFine |
number of discretization steps for the fine trajectory. Must
be larger than |
deltaFine |
discretization step for the fine trajectory. Must be
smaller than |
The fine trajectory is subsampled using the indexes
seq(1, NFine + 1, by = NFine / N).
A matrix of size c(N + 1, 2) containing x0 in the
first entry and the discretized trajectory.
samp <- rTrajWn2D(x0 = c(0, 0), alpha = c(1, 1, -0.5), mu = c(pi, pi), sigma = c(1, 1), N = 1000, delta = 0.01) plot(samp, xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 19, cex = 0.25, xlab = expression(X[t]), ylab = expression(Y[t]), col = rainbow(1000)) linesTorus(samp[, 1], samp[, 2], col = rainbow(1000))samp <- rTrajWn2D(x0 = c(0, 0), alpha = c(1, 1, -0.5), mu = c(pi, pi), sigma = c(1, 1), N = 1000, delta = 0.01) plot(samp, xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 19, cex = 0.25, xlab = expression(X[t]), ylab = expression(Y[t]), col = rainbow(1000)) linesTorus(samp[, 1], samp[, 2], col = rainbow(1000))
Computes the weights from ,
in a safe way to avoid overflows and to truncate automatically to zero low values of .
safeSoftMax(logs, expTrc = 30)safeSoftMax(logs, expTrc = 30)
logs |
matrix of logarithms where each row contains a set of |
expTrc |
truncation for exponential: |
The logs argument must be always a matrix.
A matrix of the size as logs containing the weights for each row.
# A matrix safeSoftMax(rbind(1:10, 20:11)) rbind(exp(1:10) / sum(exp(1:10)), exp(20:11) / sum(exp(20:11))) # A row-matrix safeSoftMax(rbind(-100:100), expTrc = 30) exp(-100:100) / sum(exp(-100:100))# A matrix safeSoftMax(rbind(1:10, 20:11)) rbind(exp(1:10) / sum(exp(1:10)), exp(20:11) / sum(exp(20:11))) # A row-matrix safeSoftMax(rbind(-100:100), expTrc = 30) exp(-100:100) / sum(exp(-100:100))
Given a wrapped normal density, find the parameters of a von Mises that matches it according to two characteristics: moments and scores. Score matching estimators are available for univariate and bivariate cases and moment matching only for the former.
scoreMatchWnBvm(Sigma = NULL, invSigma) scoreMatchWnVm(sigma, sigma2 = NULL) momentMatchWnVm(sigma, sigma2 = NULL)scoreMatchWnBvm(Sigma = NULL, invSigma) scoreMatchWnVm(sigma, sigma2 = NULL) momentMatchWnVm(sigma, sigma2 = NULL)
Sigma, invSigma
|
covariance or precision matrix of the bivariate wrapped normal. |
sigma, sigma2
|
standard deviation or variance of the wrapped normal. |
If the precision matrix is singular or if there are no solutions for
the score matching estimator, c(0, 0, 0) is returned.
Vector of parameters , where
is a suitable input for kappa in
dBvm.
Mardia, K. V., Kent, J. T., and Laha, A. K. (2016). Score matching estimators for directional distributions. arXiv:1604.0847. https://arxiv.org/abs/1604.08470
# Univariate WN approximation sigma <- 0.5 curve(dWn1D(x = x, mu = 0, sigma = sigma), from = -pi, to = pi, ylab = "Density", ylim = c(0, 1)) curve(dVm(x = x, mu = 0, kappa = momentMatchWnVm(sigma = sigma)), from = -pi, to = pi, col = "red", add = TRUE) curve(dVm(x = x, mu = 0, kappa = scoreMatchWnVm(sigma = sigma)), from = -pi, to = pi, col = "green", add = TRUE) # Bivariate WN approximation # WN alpha <- c(2, 1, 1) sigma <- c(1, 1) mu <- c(pi / 2, pi / 2) x <- seq(-pi, pi, l = 101)[-101] plotSurface2D(x, x, f = function(x) dStatWn2D(x = x, alpha = alpha, mu = mu, sigma = sigma), fVect = TRUE) A <- alphaToA(alpha = alpha, sigma = sigma) S <- 0.5 * solve(A) %*% diag(sigma) # Score matching kappa <- scoreMatchWnBvm(Sigma = S) # dBvm uses lambda / 2 in the exponent plotSurface2D(x, x, f = function(x) dBvm(x = x, mu = mu, kappa = c(kappa[1:2], 2 * kappa[3])), fVect = TRUE) # With singular Sigma invSigma <- matrix(c(1, sqrt(0.999), sqrt(0.999), 1), nrow = 2, ncol = 2) scoreMatchWnBvm(invSigma = invSigma) invSigma <- matrix(1, nrow = 2, ncol = 2) scoreMatchWnBvm(invSigma = invSigma)# Univariate WN approximation sigma <- 0.5 curve(dWn1D(x = x, mu = 0, sigma = sigma), from = -pi, to = pi, ylab = "Density", ylim = c(0, 1)) curve(dVm(x = x, mu = 0, kappa = momentMatchWnVm(sigma = sigma)), from = -pi, to = pi, col = "red", add = TRUE) curve(dVm(x = x, mu = 0, kappa = scoreMatchWnVm(sigma = sigma)), from = -pi, to = pi, col = "green", add = TRUE) # Bivariate WN approximation # WN alpha <- c(2, 1, 1) sigma <- c(1, 1) mu <- c(pi / 2, pi / 2) x <- seq(-pi, pi, l = 101)[-101] plotSurface2D(x, x, f = function(x) dStatWn2D(x = x, alpha = alpha, mu = mu, sigma = sigma), fVect = TRUE) A <- alphaToA(alpha = alpha, sigma = sigma) S <- 0.5 * solve(A) %*% diag(sigma) # Score matching kappa <- scoreMatchWnBvm(Sigma = S) # dBvm uses lambda / 2 in the exponent plotSurface2D(x, x, f = function(x) dBvm(x = x, mu = mu, kappa = c(kappa[1:2], 2 * kappa[3])), fVect = TRUE) # With singular Sigma invSigma <- matrix(c(1, sqrt(0.999), sqrt(0.999), 1), nrow = 2, ncol = 2) scoreMatchWnBvm(invSigma = invSigma) invSigma <- matrix(1, nrow = 2, ncol = 2) scoreMatchWnBvm(invSigma = invSigma)
Implementation of statistical methods for the estimation of toroidal diffusions. Several diffusive models are provided, most of them belonging to the Langevin family of diffusions on the torus. Specifically, the wrapped normal and von Mises processes are included, which can be seen as toroidal analogues of the Ornstein–Uhlenbeck diffusion. A collection of methods for approximate maximum likelihood estimation, organized in four blocks, is given: (i) based on the exact transition probability density, obtained as the numerical solution to the Fokker-Plank equation; (ii) based on wrapped pseudo-likelihoods; (iii) based on specific analytic approximations by wrapped processes; (iv) based on maximum likelihood of the stationary densities. The package allows the replicability of the results in García-Portugués et al. (2019) <doi:10.1007/s11222-017-9790-2>.
Eduardo García-Portugués.
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
Estimation of the in the multivariate diffusion
by the high-frequency estimate
sigmaDiff(data, delta, circular = TRUE, diagonal = FALSE, isotropic = FALSE)sigmaDiff(data, delta, circular = TRUE, diagonal = FALSE, isotropic = FALSE)
data |
vector or matrix of size |
delta |
discretization step. |
circular |
whether the process is circular or not. |
diagonal, isotropic
|
enforce different constraints for the diffusion matrix. |
See Section 3.1 in García-Portugués et al. (2019) for details.
The estimated diffusion matrix of size c(p, p).
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
# 1D x <- drop(euler1D(x0 = 0, alpha = 1, mu = 0, sigma = 1, N = 1000, delta = 0.01)) sigmaDiff(x, delta = 0.01) # 2D x <- t(euler2D(x0 = rbind(c(pi, pi)), A = rbind(c(2, 1), c(1, 2)), mu = c(pi, pi), sigma = c(1, 1), N = 1000, delta = 0.01)[1, , ]) sigmaDiff(x, delta = 0.01) sigmaDiff(x, delta = 0.01, circular = FALSE) sigmaDiff(x, delta = 0.01, diagonal = TRUE) sigmaDiff(x, delta = 0.01, isotropic = TRUE)# 1D x <- drop(euler1D(x0 = 0, alpha = 1, mu = 0, sigma = 1, N = 1000, delta = 0.01)) sigmaDiff(x, delta = 0.01) # 2D x <- t(euler2D(x0 = rbind(c(pi, pi)), A = rbind(c(2, 1), c(1, 2)), mu = c(pi, pi), sigma = c(1, 1), N = 1000, delta = 0.01)[1, , ]) sigmaDiff(x, delta = 0.01) sigmaDiff(x, delta = 0.01, circular = FALSE) sigmaDiff(x, delta = 0.01, diagonal = TRUE) sigmaDiff(x, delta = 0.01, isotropic = TRUE)
Implementation of the Thomas algorithm to solve efficiently the tridiagonal matrix system
with (usual tridiagonal matrix). If or (circulant tridiagonal matrix), then the Sherman–Morrison formula is employed.
solveTridiag(a, b, c, d, LU = 0L) solveTridiagMatConsts(a, b, c, d, LU = 0L) solvePeriodicTridiag(a, b, c, d, LU = 0L) forwardSweepTridiag(a, b, c) forwardSweepPeriodicTridiag(a, b, c)solveTridiag(a, b, c, d, LU = 0L) solveTridiagMatConsts(a, b, c, d, LU = 0L) solvePeriodicTridiag(a, b, c, d, LU = 0L) forwardSweepTridiag(a, b, c) forwardSweepPeriodicTridiag(a, b, c)
a, b, c
|
subdiagonal (below main diagonal), diagonal and superdiagonal (above main diagonal), respectively. They all are vectors of length |
d |
vector of constant terms, of length |
LU |
flag denoting if the forward sweep encoding the LU decomposition is supplied in vectors |
The Thomas algorithm is stable if the matrix is diagonally dominant.
For the periodic case, two non-periodic tridiagonal systems with different constant terms (but same coefficients) are solved using solveTridiagMatConsts. These two solutions are combined by the Sherman–Morrison formula to obtain the solution to the periodic system.
Note that the output of solveTridiag and solveTridiagMatConsts are independent from the values of a[1] and c[n], but solvePeriodicTridiag is not.
If LU is TRUE, then b and c must be precomputed with forwardSweepTridiag orforwardSweepPeriodicTridiag for its use in the call of the appropriate solver, which will be slightly faster.
solve* functions: the solution, a vector of length n and a matrix with n rows forsolveTridiagMatConsts.
forward* functions: the matrix cbind(b, c) creating the suitable b and c arguments for calling solve* when LU is TRUE.
Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Springer, New York. doi:10.1007/978-1-4899-7278-1
Conte, S. D. and de Boor, C. (1980). Elementary Numerical Analysis: An Algorithmic Approach. Third edition. McGraw-Hill, New York. doi:10.1137/1.9781611975208
# Tridiagonal matrix n <- 10 a <- rnorm(n, 3, 1) b <- rnorm(n, 10, 1) c <- rnorm(n, 0, 1) d <- rnorm(n, 0, 1) A <- matrix(0, nrow = n, ncol = n) diag(A) <- b for (i in 1:(n - 1)) { A[i + 1, i] <- a[i + 1] A[i, i + 1] <- c[i] } A # Same solutions drop(solveTridiag(a = a, b = b, c = c, d = d)) solve(a = A, b = d) # Presaving the forward sweep (encodes the LU factorization) LU <- forwardSweepTridiag(a = a, b = b, c = c) drop(solveTridiag(a = a, b = LU[, 1], c = LU[, 2], d = d, LU = 1)) # With equal coefficient matrix solveTridiagMatConsts(a = a, b = b, c = c, d = cbind(d, d + 1)) cbind(solve(a = A, b = d), solve(a = A, b = d + 1)) LU <- forwardSweepTridiag(a = a, b = b, c = c) solveTridiagMatConsts(a = a, b = LU[, 1], c = LU[, 2], d = cbind(d, d + 1), LU = 1) # Periodic matrix A[1, n] <- a[1] A[n, 1] <- c[n] A # Same solutions drop(solvePeriodicTridiag(a = a, b = b, c = c, d = d)) solve(a = A, b = d) # Presaving the forward sweep (encodes the LU factorization) LU <- forwardSweepPeriodicTridiag(a = a, b = b, c = c) drop(solvePeriodicTridiag(a = a, b = LU[, 1], c = LU[, 2], d = d, LU = 1))# Tridiagonal matrix n <- 10 a <- rnorm(n, 3, 1) b <- rnorm(n, 10, 1) c <- rnorm(n, 0, 1) d <- rnorm(n, 0, 1) A <- matrix(0, nrow = n, ncol = n) diag(A) <- b for (i in 1:(n - 1)) { A[i + 1, i] <- a[i + 1] A[i, i + 1] <- c[i] } A # Same solutions drop(solveTridiag(a = a, b = b, c = c, d = d)) solve(a = A, b = d) # Presaving the forward sweep (encodes the LU factorization) LU <- forwardSweepTridiag(a = a, b = b, c = c) drop(solveTridiag(a = a, b = LU[, 1], c = LU[, 2], d = d, LU = 1)) # With equal coefficient matrix solveTridiagMatConsts(a = a, b = b, c = c, d = cbind(d, d + 1)) cbind(solve(a = A, b = d), solve(a = A, b = d + 1)) LU <- forwardSweepTridiag(a = a, b = b, c = c) solveTridiagMatConsts(a = a, b = LU[, 1], c = LU[, 2], d = cbind(d, d + 1), LU = 1) # Periodic matrix A[1, n] <- a[1] A[n, 1] <- c[n] A # Same solutions drop(solvePeriodicTridiag(a = a, b = b, c = c, d = d)) solve(a = A, b = d) # Presaving the forward sweep (encodes the LU factorization) LU <- forwardSweepPeriodicTridiag(a = a, b = b, c = c) drop(solvePeriodicTridiag(a = a, b = LU[, 1], c = LU[, 2], d = d, LU = 1))
Simulates M trajectories starting from different initial values x0 of the WN or vM diffusion in 1D, by the Euler method, and returns their ends.
stepAheadWn1D(x0, alpha, mu, sigma, M, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)stepAheadWn1D(x0, alpha, mu, sigma, M, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)
x0 |
vector of length |
alpha |
drift parameter. |
mu |
mean parameter. Must be in |
sigma |
diffusion coefficient. |
M |
number of Monte Carlo replicates. |
N |
number of discretization steps. |
delta |
discretization step. |
type |
integer giving the type of diffusion. Currently, only |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
A matrix of size c(nx0, M) containing the M trajectory ends for each starting value x0.
N <- 100 nx0 <- 20 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] set.seed(12345678) samp1 <- euler1D(x0 = x0, mu = 0, alpha = 3, sigma = 1, N = N, delta = 0.01, type = 2) tt <- seq(0, 1, l = N + 1) plot(rep(0, nx0), x0, pch = 16, col = rainbow(nx0), xlim = c(0, 1)) for (i in 1:nx0) linesCirc(tt, samp1[i, ], col = rainbow(nx0)[i]) set.seed(12345678) samp2 <- stepAheadWn1D(x0 = x0, mu = 0, alpha = 3, sigma = 1, M = 1, N = N, delta = 0.01, type = 2) points(rep(1, nx0), samp2[, 1], pch = 16, col = rainbow(nx0)) samp1[, N + 1] samp2[, 1]N <- 100 nx0 <- 20 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] set.seed(12345678) samp1 <- euler1D(x0 = x0, mu = 0, alpha = 3, sigma = 1, N = N, delta = 0.01, type = 2) tt <- seq(0, 1, l = N + 1) plot(rep(0, nx0), x0, pch = 16, col = rainbow(nx0), xlim = c(0, 1)) for (i in 1:nx0) linesCirc(tt, samp1[i, ], col = rainbow(nx0)[i]) set.seed(12345678) samp2 <- stepAheadWn1D(x0 = x0, mu = 0, alpha = 3, sigma = 1, M = 1, N = N, delta = 0.01, type = 2) points(rep(1, nx0), samp2[, 1], pch = 16, col = rainbow(nx0)) samp1[, N + 1] samp2[, 1]
Simulates M trajectories starting from different initial values x0 of the WN or MvM diffusion in 2D, by the Euler method, and returns their ends.
stepAheadWn2D(x0, mu, A, sigma, rho = 0, M = 100L, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)stepAheadWn2D(x0, mu, A, sigma, rho = 0, M = 100L, N = 100L, delta = 0.01, type = 1L, maxK = 2L, expTrc = 30)
x0 |
matrix of size |
mu |
a vector of length |
A |
drift matrix of size |
sigma |
vector of length |
rho |
correlation coefficient of |
M |
number of Monte Carlo replicates. |
N |
number of discretization steps. |
delta |
discretization step. |
type |
integer giving the type of diffusion. Currently, only |
maxK |
maximum absolute value of the windings considered in the computation of the WN. |
expTrc |
truncation for exponential: |
An array of size c(nx0, 2, M) containing the M trajectory ends for each starting value x0.
N <- 100 nx0 <- 3 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] x0 <- as.matrix(expand.grid(x0, x0)) nx0 <- nx0^2 set.seed(12345678) samp1 <- euler2D(x0 = x0, mu = c(0, 0), A = rbind(c(3, 1), 1:2), sigma = c(1, 1), N = N, delta = 0.01, type = 2) plot(x0[, 1], x0[, 2], xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 16, col = rainbow(nx0)) for (i in 1:nx0) linesTorus(samp1[i, 1, ], samp1[i, 2, ], col = rainbow(nx0, alpha = 0.75)[i]) set.seed(12345678) samp2 <- stepAheadWn2D(x0 = x0, mu = c(0, 0), A = rbind(c(3, 1), 1:2), sigma = c(1, 1), M = 2, N = N, delta = 0.01, type = 2) points(samp2[, 1, 1], samp2[, 2, 1], pch = 16, col = rainbow(nx0)) samp1[, , N + 1] samp2[, , 1]N <- 100 nx0 <- 3 x0 <- seq(-pi, pi, l = nx0 + 1)[-(nx0 + 1)] x0 <- as.matrix(expand.grid(x0, x0)) nx0 <- nx0^2 set.seed(12345678) samp1 <- euler2D(x0 = x0, mu = c(0, 0), A = rbind(c(3, 1), 1:2), sigma = c(1, 1), N = N, delta = 0.01, type = 2) plot(x0[, 1], x0[, 2], xlim = c(-pi, pi), ylim = c(-pi, pi), pch = 16, col = rainbow(nx0)) for (i in 1:nx0) linesTorus(samp1[i, 1, ], samp1[i, 2, ], col = rainbow(nx0, alpha = 0.75)[i]) set.seed(12345678) samp2 <- stepAheadWn2D(x0 = x0, mu = c(0, 0), A = rbind(c(3, 1), 1:2), sigma = c(1, 1), M = 2, N = N, delta = 0.01, type = 2) points(samp2[, 1, 1], samp2[, 2, 1], pch = 16, col = rainbow(nx0)) samp1[, , N + 1] samp2[, , 1]
Utilities for transforming a reals into ,
or .
toPiInt(x) to2PiInt(x) toInt(x, a, b)toPiInt(x) to2PiInt(x) toInt(x, a, b)
x |
a vector, matrix or object for whom |
a, b
|
the lower and upper limits of |
Note that is excluded from the result, see examples.
The wrapped vector in the chosen interval.
# Wrapping of angles x <- seq(-3 * pi, 5 * pi, l = 100) toPiInt(x) to2PiInt(x) # Transformation to [1, 5) x <- 1:10 toInt(x, 1, 5) toInt(x + 1, 1, 5) # Transformation to [1, 5] toInt(x, 1, 6) toInt(x + 1, 1, 6)# Wrapping of angles x <- seq(-3 * pi, 5 * pi, l = 100) toPiInt(x) to2PiInt(x) # Transformation to [1, 5) x <- 1:10 toInt(x, 1, 5) toInt(x + 1, 1, 5) # Transformation to [1, 5] toInt(x, 1, 6) toInt(x + 1, 1, 6)
Wrapper for drawing pretty axis labels for circular variables.
To be invoked after plot with axes = FALSE has been called.
torusAxis(sides = 1:2, twoPi = FALSE, ...)torusAxis(sides = 1:2, twoPi = FALSE, ...)
sides |
an integer vector specifying which side of the plot the axes are
to be drawn on. The axes are placed as follows: |
twoPi |
flag indicating that |
... |
further parameters passed to |
The function calls box.
This function is usually invoked for its side effect, which is to add axes to an already existing plot.
grid <- seq(-pi, pi, l = 100) plotSurface2D(grid, grid, f = function(x) sin(x[1]) * cos(x[2]), nLev = 20, axes = FALSE) torusAxis()grid <- seq(-pi, pi, l = 100) plotSurface2D(grid, grid, f = function(x) sin(x[1]) * cos(x[2]), nLev = 20, axes = FALSE) torusAxis()
Wrapper for drawing pretty axis labels for circular variables.
To be invoked after plot3d with axes = FALSE and
box = FALSE has been called.
torusAxis3d(sides = 1:3, twoPi = FALSE, ...)torusAxis3d(sides = 1:3, twoPi = FALSE, ...)
sides |
an integer vector specifying which side of the plot the axes are
to be drawn on. The axes are placed as follows: |
twoPi |
flag indicating that |
... |
further parameters passed to |
The function calls box3d.
This function is usually invoked for its side effect, which is to add axes to an already existing plot.
if (requireNamespace("rgl")) { n <- 50 x <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) y <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) z <- toPiInt(x + y + rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) rgl::plot3d(x, y, z, xlim = c(-pi, pi), ylim = c(-pi, pi), zlim = c(-pi, pi), col = rainbow(n), size = 2, box = FALSE, axes = FALSE) torusAxis3d() }if (requireNamespace("rgl")) { n <- 50 x <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) y <- toPiInt(rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) z <- toPiInt(x + y + rnorm(n, mean = seq(-pi, pi, l = n), sd = 0.5)) rgl::plot3d(x, y, z, xlim = c(-pi, pi), ylim = c(-pi, pi), zlim = c(-pi, pi), col = rainbow(n), size = 2, box = FALSE, axes = FALSE) torusAxis3d() }
Completes a circular time series to a linear one by computing the closest wind numbers. Useful for plotting circular trajectories with crossing of boundaries.
unwrapCircSeries(x)unwrapCircSeries(x)
x |
wrapped angles. Must be a vector or a matrix, see details. |
If x is a matrix then the unwrapping is carried out
row-wise, on each column separately.
A vector or matrix containing a choice of unwrapped angles of
x that maximizes linear continuity.
# Vectors x <- c(-pi, -pi/2, pi - 0.1, -pi + 0.2) u <- unwrapCircSeries(x) max(abs(toPiInt(u) - x)) plot(toPiInt(x), ylim = c(-pi, pi)) for(k in -1:1) lines(u + 2 * k * pi) # Matrices N <- 100 set.seed(234567) x <- t(euler2D(x0 = rbind(c(0, 0)), A = diag(c(1, 1)), sigma = rep(1, 2), mu = c(pi, pi), N = N, delta = 1, type = 2)[1, , ]) u <- unwrapCircSeries(x) max(abs(toPiInt(u) - x)) plot(x, xlim = c(-pi, pi), ylim = c(-pi, pi)) for(k1 in -3:3) for(k2 in -3:3) lines(u[, 1] + 2 * k1 * pi, u[, 2] + 2 * k2 * pi, col = gray(0.5))# Vectors x <- c(-pi, -pi/2, pi - 0.1, -pi + 0.2) u <- unwrapCircSeries(x) max(abs(toPiInt(u) - x)) plot(toPiInt(x), ylim = c(-pi, pi)) for(k in -1:1) lines(u + 2 * k * pi) # Matrices N <- 100 set.seed(234567) x <- t(euler2D(x0 = rbind(c(0, 0)), A = diag(c(1, 1)), sigma = rep(1, 2), mu = c(pi, pi), N = N, delta = 1, type = 2)[1, , ]) u <- unwrapCircSeries(x) max(abs(toPiInt(u) - x)) plot(x, xlim = c(-pi, pi), ylim = c(-pi, pi)) for(k1 in -3:3) for(k2 in -3:3) lines(u[, 1] + 2 * k1 * pi, u[, 2] + 2 * k2 * pi, col = gray(0.5))