DE.FEM.time {fdaPDE} | R Documentation |
Nonparametric spatio-temporal density estimation with differential regularization
Description
This function implements a nonparametric spatio-temporal density estimation method with differential regularization (given by the sum of the square of the L2 norm of the laplacian of the density function and the square of the L2 norm of the second- order time-derivative), when points are located over a planar mesh. The computation relies only on the C++ implementation of the algorithm.
Usage
DE.FEM.time(data, data_time, FEMbasis, mesh_time, lambda, lambda_time, scaling=NULL,
fvec=NULL, heatStep=0.1, heatIter=10, stepProposals=NULL, tol1=1e-4,
tol2=0, print=FALSE, nfolds=NULL, nsimulations=500,
step_method="Fixed_Step", direction_method="BFGS",
preprocess_method="NoCrossValidation", search="tree",
isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE,
inference = FALSE)
Arguments
data |
A matrix of dimensions #observations-by-ndim. Data are locations: each row corresponds to one point,
the first column corresponds to the |
data_time |
A vector of length #observations. The i-th datum is the time instant during which the i-th location is observed (according to the order in which locations are provided in data). |
FEMbasis |
A |
mesh_time |
A vector containing the b-splines knots for separable smoothing. It is the time mesh of the considered time domain (interval). Its nodes are in increasing order. |
lambda |
A scalar or vector of smoothing parameters in space. If it is a vector, the optimal smoothing parameter in space
is chosen, together with the optimal smoothing parameter in time, with a |
lambda_time |
A scalar or vector of smoothing parameters in time. If it is a vector, the optimal smoothing parameter in time
is chosen, together with the optimal smoothing parameter in space, with a |
scaling |
A positive factor needed to scale the smoothing parameters in the construction of confidence intervals. If the scaling is not specified, it is automatically set as the square root of the number of observations. |
fvec |
A vector of length # |
heatStep |
A real specifying the time step for the discretized heat diffusion process. Default is |
heatIter |
An integer specifying the number of iterations to perform the discretized heat diffusion process. Default is |
stepProposals |
A scalar or a vector containing the step parameters useful for the descent algorithm. If there is a
vector of parameters, the biggest one such that the functional decreases at each iteration is chosen. If it is |
tol1 |
A scalar specifying the tolerance to use for the termination criterion based on the percentage difference between two consecutive iterations of the minimization algorithm of the loss function, the log-likelihood and the penalizations. Default is 1e-5. |
tol2 |
A scalar specifying the tolerance to use for the termination criterion based on the norm of the gradient of the functional to be minimized (the true minimum is such that this norm is zero). The default version does not use this criterion. Default is 0. |
print |
A boolean that is |
nfolds |
An integer specifying the number of folds used in cross validation technique to find the best pair of
( |
nsimulations |
An integer specifying the number of iterations used in the optimization algorithms. Default value is 500. |
step_method |
A string specifying which step method to use in the descent algorithm.
If it is |
direction_method |
A string specifying which descent direction to use in the descent algorithm.
If it is |
preprocess_method |
A string specifying the k fold cross validation technique to use, if there is more than one pair of
smoothing parameters in space and in time ( |
search |
A flag to decide the search algorithm type (tree or naive or walking search algorithm). Default is |
isTimeDiscrete |
A boolean specifying the time data type: |
flagMass |
A boolean specifying whether to consider full mass matrices ( |
flagLumped |
A boolean specifying whether to perform mass lumping. This numerical technique presents computational
advantages during the procedure involving a mass matrix inversion for the computation of the space penalty matrix.
Default is |
inference |
A boolean that is |
Value
A list with the following variables:
FEMbasis |
Given FEMbasis with tree information. |
g |
A vector of length # |
f_init |
A # |
lambda |
A scalar representing the optimal smoothing parameter in space selected, together with |
lambda_time |
A scalar representing the optimal smoothing parameter in time selected, together with |
data |
A matrix of dimensions #observations-by-ndim containing the spatial data used in the algorithm. They are the same given in input if the domain is 2D pr 3D; they are the original data projected on the mesh if the domain is 2.5D. Data lying outside the spatial domain, defined through its mesh, are not considered. |
data_time |
A vector of length #observations containing the time data used in the algorithm. Data lying outside the temporal domain, defined through its mesh, are not considered. |
CV_err |
A vector of length |
Examples
library(fdaPDE)
## Create a 2D mesh over a squared domain
Xbound <- seq(-3, 3, length.out = 10)
Ybound <- seq(-3, 3, length.out = 10)
grid_XY <- expand.grid(Xbound, Ybound)
Bounds <- grid_XY[(grid_XY$Var1 %in% c(-3, 3)) | (grid_XY$Var2 %in% c(-3, 3)), ]
mesh <- create.mesh.2D(nodes = Bounds, order = 1)
mesh <- refine.mesh.2D(mesh, maximum_area = 0.25, minimum_angle = 20)
FEMbasis <- create.FEM.basis(mesh)
## Create a 1D time mesh over a (non-negative) interval
mesh_time <- seq(0, 1, length.out=3)
## Generate data
n <- 50
set.seed(10)
x <- rnorm(n,0,2)
y <- rnorm(n,0,2)
locations <- cbind(x,y)
times <- runif(n,0,1)
data <- cbind(locations, times)
t <- 0.5 # time instant in which to evaluate the solution
plot(mesh)
sample <- data[abs(data[,3]-t)<0.05,1:2]
points(sample, col="red", pch=19, cex=1, main=paste('Sample | ', t-0.05,' < t < ', t+0.05))
## Density Estimation
lambda <- 0.1
lambda_time <- 0.001
sol <- DE.FEM.time(data = locations, data_time = times, FEMbasis = FEMbasis, mesh_time = mesh_time,
lambda = lambda, lambda_time = lambda_time, nsimulations=300, inference=TRUE)
## Visualization
n = 100
X <- seq(-3, 3, length.out = n)
Y <- seq(-3, 3, length.out = n)
grid <- expand.grid(X, Y)
FEMfunction = FEM.time(sol$g, mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation <- eval.FEM.time(FEM.time = FEMfunction, locations = grid, time.instants = t)
FEMfunction_L = FEM.time(sol$g_CI_L, mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation_L <- eval.FEM.time(FEM.time = FEMfunction_L, locations = grid, time.instants = t)
FEMfunction_U = FEM.time(sol$g_CI_U, mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation_U <- eval.FEM.time(FEM.time = FEMfunction_U, locations = grid, time.instants = t)
image2D(x = X, y = Y, z = matrix(exp(evaluation_L), n, n), col = heat.colors(100),
xlab = "x", ylab = "y", contour = list(drawlabels = FALSE),
main = paste("Estimated CI lower bound at t = ", t), zlim=c(0,0.3), asp = 1)
image2D(x = X, y = Y, z = matrix(exp(evaluation), n, n), col = heat.colors(100),
xlab = "x", ylab = "y", contour = list(drawlabels = FALSE),
main = paste("Estimated density at t = ", t), zlim=c(0,0.3), asp = 1)
image2D(x = X, y = Y, z = matrix(exp(evaluation_U), n, n), col = heat.colors(100),
xlab = "x", ylab = "y", contour = list(drawlabels = FALSE),
main = paste("Estimated CI upper bound at t = ", t), zlim=c(0,0.3), asp = 1)