| Title: | Physics-Informed Spatial and Functional Data Analysis |
|---|---|
| Description: | An implementation of regression models with partial differential regularizations, making use of the Finite Element Method. The models efficiently handle data distributed over irregularly shaped domains and can comply with various conditions at the boundaries of the domain. A priori information about the spatial structure of the phenomenon under study can be incorporated in the model via the differential regularization. See Sangalli, L. M. (2021) <doi:10.1111/insr.12444> "Spatial Regression With Partial Differential Equation Regularisation" for an overview. The release 1.1-9 requires R (>= 4.2.0) to be installed on windows machines. |
| Authors: | Eleonora Arnone [aut, cre], Aldo Clemente [aut], Laura M. Sangalli [aut], Eardi Lila [aut], Jim Ramsay [aut], Luca Formaggia [aut], Giovanni Ardenghi [ctb], Blerta Begu [ctb], Michele Cavazzutti [ctb], Alessandra Colli [ctb], Alberto Colombo [ctb], Luca Colombo [ctb], Carlo de Falco [ctb], Enrico Dall'Acqua [ctb], Giulia Ferla [ctb], Lorenzo Ghilotti [ctb], Cristina Galimberti [ctb], Jiyoung Kim [ctb], Martina Massardi [ctb], Giorgio Meretti [ctb], Simone Panzeri [ctb], Giulio Perin [ctb], Clara Pigolotti [ctb], Andrea Poiatti [ctb], Gian Matteo Rinaldi [ctb], Stefano Spaziani [ctb], Andrea Vicini [ctb], David C. Sterratt [cph] (Author of included Triangle source files) |
| Maintainer: | Eleonora Arnone <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1-21 |
| Built: | 2025-03-09 05:42:11 UTC |
| Source: | https://github.com/cran/fdaPDE |
Implements a finite area test function the horseshoe domain.
covs.test(x, y)covs.test(x, y)
x, y
|
Points at which to evaluate the test function. |
Returns function evaluations.
Sets up a Finite Element basis. It requires a mesh.2D, mesh.2.5D or mesh.3D object,
as input.
The basis' functions are globally continuos functions, that are polynomials once restricted to a triangle in the mesh.
The current implementation includes linear finite elements (when order = 1 in the input mesh) and
quadratic finite elements (when order = 2 in the input mesh).
If saveTree flag is TRUE, it saves the tree mesh information in advance inside mesh object and can be used later on to save mesh construction time.
create.FEM.basis(mesh, saveTree = FALSE)create.FEM.basis(mesh, saveTree = FALSE)
mesh |
A |
saveTree |
a flag to decide to save the tree mesh information in advance (default is FALSE) |
A FEMbasis object. This contains the mesh, along with some additional quantities:
orderEither "1" or "2" for the 2D and 2.5D case, and "1" for the 3D case. Order of the Finite Element basis.
nbasisScalar. The number of basis.
create.mesh.2D, create.mesh.2.5D,create.mesh.3D
library(fdaPDE) ## Upload the quasicircle2D data data(quasicircle2D) ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(quasicircle2D$boundary_nodes, quasicircle2D$locations), segments = quasicircle2D$boundary_segments) ## Plot it plot(mesh) ## Create the basis FEMbasis = create.FEM.basis(mesh) ## Upload the hub2.5D data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create the 2.5D mesh mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) ## Plot it plot(mesh) ## Create the basis FEMbasis = create.FEM.basis(mesh)library(fdaPDE) ## Upload the quasicircle2D data data(quasicircle2D) ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(quasicircle2D$boundary_nodes, quasicircle2D$locations), segments = quasicircle2D$boundary_segments) ## Plot it plot(mesh) ## Create the basis FEMbasis = create.FEM.basis(mesh) ## Upload the hub2.5D data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create the 2.5D mesh mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) ## Plot it plot(mesh) ## Create the basis FEMbasis = create.FEM.basis(mesh)
Create a 1.5D linear network mesh
create.mesh.1.5D(nodes, edges = NULL, order = 1, nodesattributes = NULL)create.mesh.1.5D(nodes, edges = NULL, order = 1, nodesattributes = NULL)
nodes |
A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes. |
edges |
A #edges-by-2 (when |
order |
Either '1' or '2'. It specifies wether each mesh should be represented by 2 nodes (the edges vertices) or by 3 nodes (the edges's vertices and midpoint).
These are
respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
nodesattributes |
A matrix with #nodes rows containing nodes' attributes. These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value of a Dirichlet boundary condition at boundary nodes added during the triangulation process. |
An object of the class mesh.1.5D with the following output:
nodesA #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
nodesmarkersA vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.
nodesattributesA matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.
edgesA #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in nodes, indicating the nodes where the edge starts from and ends to.
neighborsA #edges-by-2 matrix of list. Each row contains the indices of the neighbouring edges. An empty entry indicates that one node of the edge is a boundary node.
orderEither '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints). These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.
mesh.2.5D object from the nodes locations and the connectivity matrixCreate a mesh.2.5D object from the nodes locations and the connectivity matrix
create.mesh.2.5D( nodes, triangles = NULL, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL )create.mesh.2.5D( nodes, triangles = NULL, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL )
nodes |
A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes. |
triangles |
A #triangles-by-3 (when |
order |
Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
These are
respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
nodesattributes |
A matrix with #nodes rows containing nodes' attributes.
These are passed unchanged to the output. This has been added for consistency with the function |
segments |
A #segments-by-2 matrix. Each row contains the row's indices in |
holes |
A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
in the triangulation, when the domain has holes. This has been added for consistency with the function |
An object of the class mesh.2.5D with the following output:
nodesA #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.
nodesmarkersA vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.
nodesattributesA matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.
trianglesA #triangles-by-3 (when order = 1) or #triangles-by-6 (when order = 2) matrix.
It specifies the triangles giving the indices in nodes of the triangles' vertices and (when nodes = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
at
https://www.cs.cmu.edu/~quake/triangle.highorder.html.
segmentsmarkersA vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in segments is a boundary segment;
an entry '0' indicates that the corresponding segment is not a boundary segment.
edgesA #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in nodes, indicating the nodes where the edge starts from and ends to.
edgesmarkersA vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in edge is a boundary edge;
an entry '0' indicates that the corresponding edge is not a boundary edge.
neighborsA #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that one edge of the triangle is a boundary edge.
holesA #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes. These are passed unchanged from the input.
orderEither '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints). These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.
library(fdaPDE) ## Upload the hub2.5D the data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create mesh from nodes and connectivity matrix: mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) plot(mesh)library(fdaPDE) ## Upload the hub2.5D the data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create mesh from nodes and connectivity matrix: mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) plot(mesh)
This function is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used
to create a triangulation of the domain of interest starting from a list of points, to be used as triangles' vertices, and a list of segments, that define the domain boundary. The resulting
mesh is a Constrained Delaunay triangulation. This is constructed in a way to preserve segments provided in the input segments without splitting them. This imput can be used to define the boundaries
of the domain. If this imput is NULL, it generates a triangulation over the
convex hull of the points.
It is also possible to create a mesh.2D from the nodes locations and the connectivity matrix.
create.mesh.2D(nodes, nodesattributes = NA, segments = NA, holes = NA, triangles = NA, order = 1, verbosity = 0)create.mesh.2D(nodes, nodesattributes = NA, segments = NA, holes = NA, triangles = NA, order = 1, verbosity = 0)
nodes |
A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes. |
nodesattributes |
A matrix with #nodes rows containing nodes' attributes. These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value of a Dirichlet boundary condition at boundary nodes added during the triangulation process. |
segments |
A #segments-by-2 matrix. Each row contains the row's indices in |
holes |
A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes. |
triangles |
A #triangles-by-3 (when |
order |
Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
These are
respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
verbosity |
This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process. When |
An object of the class mesh.2D with the following output:
nodesA #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
nodesmarkersA vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.
nodesattributesA matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.
trianglesA #triangles-by-3 (when order = 1) or #triangles-by-6 (when order = 2) matrix.
This option is used when a triangulation is already available. It specifies the triangles giving the indices in nodes of the triangles' vertices and (when nodes = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
at
https://www.cs.cmu.edu/~quake/triangle.highorder.html.
segmentsmarkersA vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in segments is a boundary segment;
an entry '0' indicates that the corresponding segment is not a boundary segment.
edgesA #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in nodes, indicating the nodes where the edge starts from and ends to.
edgesmarkersA vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in edge is a boundary edge;
an entry '0' indicates that the corresponding edge is not a boundary edge.
neighborsA #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that one edge of the triangle is a boundary edge.
holesA #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes.
orderEither '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints). These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.
refine.mesh.2D, create.FEM.basis
library(fdaPDE) ## Upload the quasicirle2D data data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data ## Create mesh from boundary ## if the domain is convex it is sufficient to call: mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations)) plot(mesh) ## if the domain is not convex, pass in addition the segments the compose the boundary: mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create mesh from data locations (without knowing the boundary) mesh = create.mesh.2D(nodes = locations) plot(mesh) ## In this case the domain is the convex hull of the data locations. ## Do this only if you do not have any information about the shape of the domain of interest.library(fdaPDE) ## Upload the quasicirle2D data data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data ## Create mesh from boundary ## if the domain is convex it is sufficient to call: mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations)) plot(mesh) ## if the domain is not convex, pass in addition the segments the compose the boundary: mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create mesh from data locations (without knowing the boundary) mesh = create.mesh.2D(nodes = locations) plot(mesh) ## In this case the domain is the convex hull of the data locations. ## Do this only if you do not have any information about the shape of the domain of interest.
mesh.3D object from the connectivity matrix and nodes locationsCreate a mesh.3D object from the connectivity matrix and nodes locations
create.mesh.3D( nodes, tetrahedrons, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL )create.mesh.3D( nodes, tetrahedrons, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL )
nodes |
A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes. |
tetrahedrons |
A #tetrahedrons-by-4 (when |
order |
Either '1' or '2'. It specifies wether each mesh tetrahedron should be represented by 4 nodes (the tetrahedron's vertices) or by 10 nodes (the tetrahedron's vertices and edge midpoints).
These are
respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
nodesattributes |
A matrix with #nodes rows containing nodes' attributes.
These are passed unchanged to the output. This has been added for consistency with the function |
segments |
A #segments-by-2 matrix. Each row contains the row's indices in |
holes |
A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
in the triangulation, when the domain has holes. This has been added for consistency with the function |
An object of the class mesh.3D with the following output:
nodesA #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.
nodesmarkersA vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.
nodesattributesA matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.
tetrahedronsA #tetrahedrons-by-4 (when order = 1) or #tetrahedrons-by-10 (when order = 2) matrix.
It specifies the tetrahedrons giving the indices in nodes of the tetrahedrons' vertices and (when nodes = 2) also if the tetrahedrons' edges midpoints.
segmentsmarkersA vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in segments is a boundary segment;
an entry '0' indicates that the corresponding segment is not a boundary segment.
facesA #faces-by-3 matrix containing all the faces of the tetrahedrons in the output triangulation. Each row contains the row's indices in nodes, indicating the nodes where the face starts from and ends to.
facesmarkersA vector of lenght #faces with entries either '1' or '0'. An entry '1' indicates that the corresponding element in faces is a boundary face;
an entry '0' indicates that the corresponding edge is not a boundary face.
neighborsA #triangles-by-4 matrix. Each row contains the indices of the four neighbouring tetrahedrons An entry '-1' indicates that one face of the tetrahedrons is a boundary face.
holesA #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes. These are passed unchanged from the input.
orderEither '1' or '2'. It specifies wether each mesh tetrahedron should be represented by 3 nodes (the tetrahedron's vertices) or by 6 nodes (the tetrahedron's vertices and midpoints). These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.
library(fdaPDE) ##Load the matrix nodes and tetrahedrons data(sphere3Ddata) nodes=sphere3Ddata$nodes tetrahedrons=sphere3Ddata$tetrahedrons ##Create the triangulated mesh from the connectivity matrix and nodes locations mesh=create.mesh.3D(nodes,tetrahedrons)library(fdaPDE) ##Load the matrix nodes and tetrahedrons data(sphere3Ddata) nodes=sphere3Ddata$nodes tetrahedrons=sphere3Ddata$tetrahedrons ##Create the triangulated mesh from the connectivity matrix and nodes locations mesh=create.mesh.3D(nodes,tetrahedrons)
This function implements a nonparametric density estimation method with differential regularization (given by the square root of the L2 norm of the laplacian of the density function), when points are located over a planar mesh. The computation relies only on the C++ implementation of the algorithm.
DE.FEM(data, FEMbasis, lambda, scaling=NULL, fvec=NULL, heatStep=0.1, heatIter=500, stepProposals=NULL,tol1=1e-4, tol2=0, print=FALSE, nfolds=0, nsimulations=500, step_method="Fixed_Step", direction_method="BFGS", preprocess_method="NoCrossValidation", search = "tree", inference = FALSE)DE.FEM(data, FEMbasis, lambda, scaling=NULL, fvec=NULL, heatStep=0.1, heatIter=500, stepProposals=NULL,tol1=1e-4, tol2=0, print=FALSE, nfolds=0, nsimulations=500, step_method="Fixed_Step", direction_method="BFGS", preprocess_method="NoCrossValidation", search = "tree", inference = FALSE)
data |
A matrix of dimensions #observations-by-ndim. Data are locations: each row corresponds to one point,
the first column corresponds to the |
FEMbasis |
A |
lambda |
A scalar or vector of smoothing parameters. If it is a vector, the optimal smoothing parameter is chosen
with a |
scaling |
A positive factor needed to scale the smoothing parameter 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 penalization. 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 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 |
nsimulations |
An integer specifying the number of iterations used in the optimization algorithms. Default value is 500. |
step_method |
String. This parameter specifies which step method use in the descent algorithm.
If it is |
direction_method |
String. This parameter specifies which descent direction use in the descent algorithm.
If it is |
preprocess_method |
String. This parameter specifies the k fold cross validation technique to use, if there is more
than one smoothing parameter |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). Default is |
inference |
A boolean that is |
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 selected via k fold cross validation, if in the input there is a vector of parameters; the scalar given in input otherwise. |
data |
A matrix of dimensions #observations-by-ndim containing the 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. |
CV_err |
A vector of length |
Ferraccioli, F., Arnone, E., Finos, L., Ramsay, J. O., Sangalli, L. M. (2021). Nonparametric density estimation over complicated domains. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 83(2), 346-368.
Arnone, E., Ferraccioli, F., Pigolotti, C., Sangalli, L.M. (2021), A roughness penalty approach to estimate densities over two-dimensional manifolds, Computational Statistics and Data Analysis, to appear.
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.2) FEMbasis <- create.FEM.basis(mesh) ## Generate data n <- 50 set.seed(10) data_x <- rnorm(n) data_y <- rnorm(n) data <- cbind(data_x, data_y) plot(mesh) points(data, col="red", pch=19, cex=0.5) ## Density Estimation lambda = 0.1 sol <- DE.FEM(data = data, FEMbasis = FEMbasis, lambda = lambda, fvec=NULL, heatStep=0.1, heatIter=500, nsimulations=300,step_method = "Fixed_Step", inference = TRUE) ## Visualization n = 100 X <- seq(-3, 3, length.out = n) Y<- seq(-3, 3, length.out = n) grid <- expand.grid(X, Y) evaluation <- eval.FEM(FEM(FEMbasis, coeff = sol$g), locations = grid) lower_bound_g <- eval.FEM(FEM(FEMbasis, coeff = sol$g_CI_L), locations = grid) upper_bound_g <- eval.FEM(FEM(FEMbasis, coeff = sol$g_CI_U), locations = grid) evaluation <- exp(evaluation) lower_bound_g <- exp(lower_bound_g) upper_bound_g <- exp(upper_bound_g) eval <- matrix(evaluation, n, n) eval_L <- matrix(lower_bound_g, n, n) eval_U <- matrix(upper_bound_g, n, n) image2D(x = X, y = Y, z = eval_L, col = heat.colors(100), xlab = "x", ylab = "y", contour = list(drawlabels = FALSE), main = "Estimated CI lower bound") image2D(x = X, y = Y, z = eval, col = heat.colors(100), xlab = "x", ylab = "y", contour = list(drawlabels = FALSE), main = "Estimated density") image2D(x = X, y = Y, z = eval_U, col = heat.colors(100), xlab = "x", ylab = "y", contour = list(drawlabels = FALSE), main = "Estimated CI upper bound")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.2) FEMbasis <- create.FEM.basis(mesh) ## Generate data n <- 50 set.seed(10) data_x <- rnorm(n) data_y <- rnorm(n) data <- cbind(data_x, data_y) plot(mesh) points(data, col="red", pch=19, cex=0.5) ## Density Estimation lambda = 0.1 sol <- DE.FEM(data = data, FEMbasis = FEMbasis, lambda = lambda, fvec=NULL, heatStep=0.1, heatIter=500, nsimulations=300,step_method = "Fixed_Step", inference = TRUE) ## Visualization n = 100 X <- seq(-3, 3, length.out = n) Y<- seq(-3, 3, length.out = n) grid <- expand.grid(X, Y) evaluation <- eval.FEM(FEM(FEMbasis, coeff = sol$g), locations = grid) lower_bound_g <- eval.FEM(FEM(FEMbasis, coeff = sol$g_CI_L), locations = grid) upper_bound_g <- eval.FEM(FEM(FEMbasis, coeff = sol$g_CI_U), locations = grid) evaluation <- exp(evaluation) lower_bound_g <- exp(lower_bound_g) upper_bound_g <- exp(upper_bound_g) eval <- matrix(evaluation, n, n) eval_L <- matrix(lower_bound_g, n, n) eval_U <- matrix(upper_bound_g, n, n) image2D(x = X, y = Y, z = eval_L, col = heat.colors(100), xlab = "x", ylab = "y", contour = list(drawlabels = FALSE), main = "Estimated CI lower bound") image2D(x = X, y = Y, z = eval, col = heat.colors(100), xlab = "x", ylab = "y", contour = list(drawlabels = FALSE), main = "Estimated density") image2D(x = X, y = Y, z = eval_U, col = heat.colors(100), xlab = "x", ylab = "y", contour = list(drawlabels = FALSE), main = "Estimated CI upper bound")
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.
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=0, nsimulations=500, step_method="Fixed_Step", direction_method="BFGS", preprocess_method="NoCrossValidation", search="tree", isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE, inference = FALSE)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=0, nsimulations=500, step_method="Fixed_Step", direction_method="BFGS", preprocess_method="NoCrossValidation", search="tree", isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE, inference = FALSE)
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 |
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 |
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)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)
This function implements two methods for the density initialization procedure.
DE.heat.FEM(data, FEMbasis, lambda=NULL, heatStep=0.1, heatIter=500, init="Heat", nFolds=5, search = "tree")DE.heat.FEM(data, FEMbasis, lambda=NULL, heatStep=0.1, heatIter=500, init="Heat", nFolds=5, search = "tree")
data |
A matrix of dimensions #observations-by-ndim. Data are locations: each row corresponds to one point,
the first column corresponds to the |
FEMbasis |
A |
lambda |
A scalar or vector of smoothing parameters. Default is NULL. It is useful only if |
heatStep |
Real specifying the time step for the discretized heat diffusionn process. |
heatIter |
Integer specifying the number of iteriations to perform the discretized heat diffusion process. |
init |
String. This parameter specifies the initialization procedure. It can be either 'Heat' or 'CV'. |
nFolds |
An integer specifying the number of folds used in cross validation techinque. It is useful only
for the case |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). |
If init = 'Heat' it returns a matrix in which each column contains the initial vector
for each lambda. If init = 'CV' it returns the initial vector associated to the lambda given.
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.2) FEMbasis <- create.FEM.basis(mesh) ## Generate data n <- 50 set.seed(10) data_x <- rnorm(n) data_y <- rnorm(n) data <- cbind(data_x, data_y) plot(mesh) points(data, col="red", pch=19, cex=0.5) ## Density initialization lambda = 0.1 sol = DE.heat.FEM(data, FEMbasis, lambda, heatStep=0.1, heatIter=500, init="Heat") ## Visualization plot(FEM(coeff=sol$f_init, FEMbasis=FEMbasis))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.2) FEMbasis <- create.FEM.basis(mesh) ## Generate data n <- 50 set.seed(10) data_x <- rnorm(n) data_y <- rnorm(n) data <- cbind(data_x, data_y) plot(mesh) points(data, col="red", pch=19, cex=0.5) ## Density initialization lambda = 0.1 sol = DE.heat.FEM(data, FEMbasis, lambda, heatStep=0.1, heatIter=500, init="Heat") ## Visualization plot(FEM(coeff=sol$f_init, FEMbasis=FEMbasis))
This function implements two methods for the density initialization procedure.
DE.heat.FEM.time(data, data_time, FEMbasis, mesh_time, lambda=NULL, lambda_time=NULL, heatStep=0.1, heatIter=10, init="Heat", nFolds=5, search="tree", isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE)DE.heat.FEM.time(data, data_time, FEMbasis, mesh_time, lambda=NULL, lambda_time=NULL, heatStep=0.1, heatIter=10, init="Heat", nFolds=5, search="tree", isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE)
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 dimensions #observations. The i-th datum is the time instant during which the i-th location is observed (according to the order in which data are provided). |
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. Default is NULL. It is useful only if |
lambda_time |
A scalar or vector of smoothing parameters in time. Default is NULL. It is useful only if |
heatStep |
A real specifying the time step for the discretized heat diffusion process. |
heatIter |
An integer specifying the number of iterations to perform the discretized heat diffusion process. |
init |
A string specifying the initialization procedure. It can be either 'Heat' or 'CV'. |
nFolds |
An integer specifying the number of folds used in the cross validation technique. It is useful only
for the case |
search |
A flag to decide the search algorithm type (tree or naive or walking search algorithm). |
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 |
If init = 'Heat' it returns a matrix in which each column contains the initial vector
for each possible pair (lambda, lambda_time). If init = 'CV' it returns the initial vector associated
to the unique pair (lambda, lambda_time) given.
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 initialization lambda = 0.1 lambda_time <- 0.001 sol = DE.heat.FEM.time(data = locations, data_time = times, FEMbasis = FEMbasis, mesh_time = mesh_time, lambda = lambda, lambda_time = lambda_time, heatStep=0.1, heatIter=10, init="Heat") ## 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$f_init[,1,1], mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE) evaluation <- eval.FEM.time(FEM.time = FEMfunction, locations = grid, time.instants = t) image2D(x = X, y = Y, z = matrix(evaluation, n, n), col = heat.colors(100), xlab = "", ylab = "", contour = list(drawlabels = FALSE), main = paste("Initial guess at t = ", t), zlim=c(0,0.2), asp = 1)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 initialization lambda = 0.1 lambda_time <- 0.001 sol = DE.heat.FEM.time(data = locations, data_time = times, FEMbasis = FEMbasis, mesh_time = mesh_time, lambda = lambda, lambda_time = lambda_time, heatStep=0.1, heatIter=10, init="Heat") ## 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$f_init[,1,1], mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE) evaluation <- eval.FEM.time(FEM.time = FEMfunction, locations = grid, time.instants = t) image2D(x = X, y = Y, z = matrix(evaluation, n, n), col = heat.colors(100), xlab = "", ylab = "", contour = list(drawlabels = FALSE), main = paste("Initial guess at t = ", t), zlim=c(0,0.2), asp = 1)
It evaluates a FEM object at the specified set of locations or areal regions. The locations are used for pointwise evaluations and incidence matrix for areal evaluations. The locations and the incidence matrix cannot be both NULL or both provided.
eval.FEM(FEM, locations = NULL, incidence_matrix = NULL, search = "tree", bary.locations = NULL)eval.FEM(FEM, locations = NULL, incidence_matrix = NULL, search = "tree", bary.locations = NULL)
FEM |
A |
locations |
A 2-columns (in 1.5D or 2D) or 3-columns (in 2.5D and 3D) matrix with the spatial locations where the FEM object should be evaluated. |
incidence_matrix |
In case of areal evaluations, the #regions-by-#elements incidence matrix defining the regions where the FEM object should be evaluated. |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). |
bary.locations |
A list with three vectors:
|
A vector or a matrix of numeric evaluations of the FEM object.
If the FEM object contains multiple finite element functions the output is a matrix, and
each row corresponds to the location (or areal region) where the evaluation has been taken, while each column
corresponds to the function evaluated.
Sangalli, L. M., Ramsay, J. O., & Ramsay, T. O. (2013). Spatial spline regression models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(4), 681-703.
Azzimonti, L., Sangalli, L. M., Secchi, P., Domanin, M., & Nobile, F. (2015). Blood flow velocity field estimation via spatial regression with PDE penalization. Journal of the American Statistical Association, 110(511), 1057-1071.
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Evaluate the finite element function in the location (1,0.5) eval.FEM(FEMfunction, locations = matrix(c(1, 0.5), ncol = 2)) ## Evaluate the mean of the finite element function over the fifth triangle of the mesh incidence_matrix = matrix(0, ncol = nrow(mesh$triangles)) incidence_matrix[1,5] = 1 eval.FEM(FEMfunction, incidence_matrix = incidence_matrix)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Evaluate the finite element function in the location (1,0.5) eval.FEM(FEMfunction, locations = matrix(c(1, 0.5), ncol = 2)) ## Evaluate the mean of the finite element function over the fifth triangle of the mesh incidence_matrix = matrix(0, ncol = nrow(mesh$triangles)) incidence_matrix[1,5] = 1 eval.FEM(FEMfunction, incidence_matrix = incidence_matrix)
It evaluates a FEM.time object at the specified set of locations or regions.
If space.time.locations is provided locations, incidence_matrix and
time.instants must be NULL. Otherwise time.instants and one of locations and
incidence_matrix must be given. In this case the evaluation is perform on the tensor grid
time.instants-by-locations (or time.instants-by-areal domains).
eval.FEM.time(FEM.time, locations = NULL, time.instants = NULL, space.time.locations = NULL, incidence_matrix = NULL, lambdaS = 1, lambdaT = 1, search = "tree", bary.locations = NULL)eval.FEM.time(FEM.time, locations = NULL, time.instants = NULL, space.time.locations = NULL, incidence_matrix = NULL, lambdaS = 1, lambdaT = 1, search = "tree", bary.locations = NULL)
FEM.time |
A |
locations |
A 2-columns (in case of planar mesh) or 3-columns(in case of 2D manifold in a 3D space or a 3D volume) matrix with the spatial locations where the FEM.time object should be evaluated. |
time.instants |
A vector with the time instants where the FEM.time object should be evaluated. |
space.time.locations |
A 3-columns (in case of planar mesh) or 4-columns(in case of 2D manifold in a 3D space or a 3D volume)
matrix with the time instants and spatial locations where the FEM.time object should be evaluated.
The first column is for the time instants. If given, |
incidence_matrix |
In case of areal data, the #regions x #elements incidence matrix defining the regions. |
lambdaS |
The index of the lambdaS choosen for the evaluation. |
lambdaT |
The index of the lambdaT choosen for the evaluation. |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). |
bary.locations |
A list with three vectors:
|
A matrix of numeric evaluations of the FEM.time object. Each row indicates the location where
the evaluation has been taken, the column indicates the function evaluated.
Devillers, O. et al. 2001. Walking in a Triangulation, Proceedings of the Seventeenth Annual Symposium on Computational Geometry
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 time = 1:5 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time ## Create the FEM.time object FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5, FEMbasis=FEMbasis, FLAG_PARABOLIC=TRUE) evaluations = eval.FEM.time(FEM_time_function, locations = matrix(c(-0.92,0), ncol=2), time.instants = time)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 time = 1:5 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time ## Create the FEM.time object FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5, FEMbasis=FEMbasis, FLAG_PARABOLIC=TRUE) evaluations = eval.FEM.time(FEM_time_function, locations = matrix(c(-0.92,0), ncol=2), time.instants = time)
Only executed when smooth.FEM.basis is run with the option CPP_CODE = FALSE. It computes the mass matrix. The element (i,j) of this matrix contains the integral over the domain of the product between the ith and kth element
of the Finite Element basis. As common practise in Finite Element Analysis, this quantities are computed iterating over all the mesh triangles.
Only executed when smooth.FEM.basis is run with the option CPP_CODE = FALSE. It computes the stifness matrix. The element (i,j) of this matrix contains the integral over the domain of the scalar product between the gradient of the ith and kth element
of the Finite Element basis. As common practise in Finite Element Analysis, this quantities are computed iterating over all the mesh triangles.
Only executed when the function smooth.FEM.basis is run with the option CPP_CODE = FALSE. It evaluates the Finite Element basis functions and their derivatives up to order 2 at the specified set of locations.
This version of the function is implemented using only R code. It is called by R_smooth.FEM.basis.
Only executed when the function smooth.FEM.basis is run with the option CPP_CODE = FALSE. It evaluates a FEM object at the specified set of locations.
This function implements a spatial regression model with differential regularization; isotropic and stationary case. In particular, the regularizing term involves the Laplacian of the spatial field. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.
This function implements a spatial regression model with differential regularization; anysotropic case. In particular, the regularizing term involves a second order elliptic PDE, that models the space-variation of the phenomenon. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.
This function implements a spatial regression model with differential regularization; anysotropic and non-stationary case. In particular, the regularizing term involves a second order elliptic PDE with space-varying coefficients, that models the space-variation of the phenomenon. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.
This function is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used
to create a triangulation of the domain of interest starting from a list of points, to be used as triangles' vertices, and a list of segments, that define the domain boundary. The resulting
mesh is a Constrained Delaunay triangulation. This is constructed in a way to preserve segments provided in the input segments without splitting them. This imput can be used to define the boundaries
of the domain. If this imput is NULL, it generates a triangulation over the
convex hull of the points.
This function refines a Constrained Delaunay triangulation into a Conforming Delaunay triangulation. This is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used to
refine a mesh created previously with create.MESH.2D. The algorithm can add Steiner points (points through which the segments are splitted)
in order to meet the imposed refinement conditions.
Plot a mesh MESH2D object, generated by create.MESH.2D or refine.MESH.2D. Circles indicate the mesh nodes.
R_mass(FEMbasis) R_stiff(FEMbasis) R_smooth.FEM.basis( locations, observations, FEMbasis, lambda, covariates = NULL, GCV ) R_eval.FEM.basis(FEMbasis, locations, nderivs = matrix(0, 1, 2)) R_eval.FEM(FEM, locations) smooth.FEM.basis( locations = NULL, observations, FEMbasis, lambda, covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE ) smooth.FEM.PDE.basis( locations = NULL, observations, FEMbasis, lambda, PDE_parameters, covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE ) smooth.FEM.PDE.sv.basis( locations = NULL, observations, FEMbasis, lambda, PDE_parameters, covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE ) create.MESH.2D(nodes, nodesattributes = NA, segments = NA, holes = NA, triangles = NA, order = 1, verbosity = 0) refine.MESH.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity) ## S3 method for class 'MESH2D' plot(x, ...)R_mass(FEMbasis) R_stiff(FEMbasis) R_smooth.FEM.basis( locations, observations, FEMbasis, lambda, covariates = NULL, GCV ) R_eval.FEM.basis(FEMbasis, locations, nderivs = matrix(0, 1, 2)) R_eval.FEM(FEM, locations) smooth.FEM.basis( locations = NULL, observations, FEMbasis, lambda, covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE ) smooth.FEM.PDE.basis( locations = NULL, observations, FEMbasis, lambda, PDE_parameters, covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE ) smooth.FEM.PDE.sv.basis( locations = NULL, observations, FEMbasis, lambda, PDE_parameters, covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE ) create.MESH.2D(nodes, nodesattributes = NA, segments = NA, holes = NA, triangles = NA, order = 1, verbosity = 0) refine.MESH.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity) ## S3 method for class 'MESH2D' plot(x, ...)
FEMbasis |
A |
locations |
A #observations-by-2 matrix where each row specifies the spatial coordinates |
observations |
A vector of length #observations with the observed data values over the domain.
The locations of the observations can be specified with the |
lambda |
A scalar or vector of smoothing parameters. |
covariates |
A #observations-by-#covariates matrix where each row represents the covariates associated with the corresponding observed data value in |
GCV |
Boolean. If |
nderivs |
A vector of lenght 2 specifying the order of the partial derivatives of the bases to be evaluated. The vectors' entries can be 0,1 or 2, where 0 indicates that only the basis functions, and not their derivatives, should be evaluated. |
FEM |
A |
BC |
A list with two vectors:
|
CPP_CODE |
Boolean. If |
PDE_parameters |
A list specifying the space-varying parameters of the elliptic PDE in the regularizing term: |
nodes |
A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes. |
nodesattributes |
A matrix with #nodes rows containing nodes' attributes. These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value of a Dirichlet boundary condition at boundary nodes added during the triangulation process. |
segments |
A #segments-by-2 matrix. Each row contains the row's indices in |
holes |
A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes. |
triangles |
A #triangles-by-3 (when |
order |
Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
These are
respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
verbosity |
This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process. |
mesh |
A MESH2D object representing the triangular mesh, created by create.MESH.2D. |
minimum_angle |
A scalar specifying a minimun value for the triangles angles. |
maximum_area |
A scalar specifying a maximum value for the triangles areas. |
delaunay |
A boolean parameter indicating whether or not the output mesh should satisfy the Delaunay condition. |
x |
A MESH2D object defining the triangular mesh, as generated by |
... |
Arguments representing graphical options to be passed to par. |
These functions are Deprecated in this release of fdaPDE, they will be marked as Defunct and removed in a future version.
A square matrix with the integrals of all the basis' functions pairwise products. The dimension of the matrix is equal to the number of the nodes of the mesh.
A square matrix with the integrals of all the basis functions' gradients pairwise dot products. The dimension of the matrix is equal to the number of the nodes of the mesh.
A list with the following quantities:
fit.FEM |
A |
PDEmisfit.FEM |
A |
beta |
If covariates is not |
edf |
If GCV is |
stderr |
If GCV is |
GCV |
If GCV is |
A matrix of basis function values. Each row indicates the location where the evaluation has been taken, the column indicates the basis function evaluated
A matrix of numeric evaluations of the FEM object. Each row indicates the location where the evaluation has been taken, the column indicates the
function evaluated.
A list with the following variables:
fit.FEMA FEM object that represents the fitted spatial field.
PDEmisfit.FEMA FEM object that represents the Laplacian of the estimated spatial field.
solutionA list, note that all terms are matrices or row vectors: the jth column represents the vector of related to lambda[j] if lambda.selection.criterion="grid" and lambda.selection.lossfunction="unused".
In all the other cases is returned just the column related to the best penalization parameter
fMatrix, estimate of function f, first half of solution vector
gMatrix, second half of solution vector
z_hatMatrix, prediction of the output in the locations
betaIf covariates is not NULL, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. It is the regression coefficients estimate
rmseEstimate of the root mean square error in the locations
estimated_sdEstiimate of the standard deviation of the error
optimizationA detailed list of optimization related data:
lambda_solutionnumerical value of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"
lambda_positioninteger, postion in lambda_vector of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"
GCVnumeric value of GCV in correspondence of the optimum
optimization_detailslist containing further information about the optimization method used and the nature of its termination, eventual number of iterations
dofnumeric vector, value of dof for all the penalizations it has been computed, empty if not computed
lambda_vectornumeric value of the penalization factors passed by the user or found in the iterations of the optimization method
GCV_vectornumeric vector, value of GCV for all the penalizations it has been computed
timeDuration of the entire optimization computation
bary.locationsA barycenter information of the given locations if the locations are not mesh nodes.
A list with the following variables:
fit.FEMA FEM object that represents the fitted spatial field.
PDEmisfit.FEMA FEM object that represents the Laplacian of the estimated spatial field.
solutionA list, note that all terms are matrices or row vectors: the jth column represents the vector of related to lambda[j] if lambda.selection.criterion="grid" and lambda.selection.lossfunction="unused".
In all the other cases is returned just the column related to the best penalization parameter
fMatrix, estimate of function f, first half of solution vector
gMatrix, second half of solution vector
z_hatMatrix, prediction of the output in the locations
betaIf covariates is not NULL, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. It is the regression coefficients estimate
rmseEstimate of the root mean square error in the locations
estimated_sdEstiimate of the standard deviation of the error
optimizationA detailed list of optimization related data:
lambda_solutionnumerical value of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"
lambda_positioninteger, postion in lambda_vector of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"
GCVnumeric value of GCV in correspondence of the optimum
optimization_detailslist containing further information about the optimization method used and the nature of its termination, eventual number of iterations
dofnumeric vector, value of dof for all the penalizations it has been computed, empty if not computed
lambda_vectornumeric value of the penalization factors passed by the user or found in the iterations of the optimization method
GCV_vectornumeric vector, value of GCV for all the penalizations it has been computed
timeDuration of the entire optimization computation
bary.locationsA barycenter information of the given locations if the locations are not mesh nodes.
A list with the following variables:
fit.FEMA FEM object that represents the fitted spatial field.
PDEmisfit.FEMA FEM object that represents the Laplacian of the estimated spatial field.
solutionA list, note that all terms are matrices or row vectors: the jth column represents the vector of related to lambda[j] if lambda.selection.criterion="grid" and lambda.selection.lossfunction="unused".
In all the other cases is returned just the column related to the best penalization parameter
fMatrix, estimate of function f, first half of solution vector
gMatrix, second half of solution vector
z_hatMatrix, prediction of the output in the locations
betaIf covariates is not NULL, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. It is the regression coefficients estimate
rmseEstimate of the root mean square error in the locations
estimated_sdEstiimate of the standard deviation of the error
optimizationA detailed list of optimization related data:
lambda_solutionnumerical value of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"
lambda_positioninteger, postion in lambda_vector of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"
GCVnumeric value of GCV in correspondence of the optimum
optimization_detailslist containing further information about the optimization method used and the nature of its termination, eventual number of iterations
dofnumeric vector, value of dof for all the penalizations it has been computed, empty if not computed
lambda_vectornumeric value of the penalization factors passed by the user or found in the iterations of the optimization method
GCV_vectornumeric vector, value of GCV for all the penalizations it has been computed
timeDuration of the entire optimization computation
bary.locationsA barycenter information of the given locations if the locations are not mesh nodes.
An object of the class MESH2D with the following output:
nodes |
A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes. |
nodesmarkers |
A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node. |
nodesattributes |
nodesattributes A matrix with #nodes rows containing nodes' attributes. These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value of a Dirichlet boundary condition at boundary nodes added during the triangulation process. |
triangles |
A #triangles-by-3 (when |
segmentsmarker |
A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in |
edges |
A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in |
edgesmarkers |
A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in |
neighbors |
A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that one edge of the triangle is a boundary edge. |
holes |
A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes. |
order |
Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
A MESH2D object representing the refined triangular mesh, with the following output:
nodes |
A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes. |
nodesmarkers |
A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node. |
nodesattributes |
nodesattributes A matrix with #nodes rows containing nodes' attributes. These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value of a Dirichlet boundary condition at boundary nodes added during the triangulation process. |
triangles |
A #triangles-by-3 (when |
edges |
A #edges-by-2 matrix. Each row contains the row's indices of the nodes where the edge starts from and ends to. |
edgesmarkers |
A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in |
neighbors |
A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that one edge of the triangle is a boundary edge. |
holes |
A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes. |
order |
Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is |
refine.MESH.2D, create.FEM.basis
create.MESH.2D, create.FEM.basis
This function defines a FEM object.
FEM(coeff,FEMbasis)FEM(coeff,FEMbasis)
coeff |
A vector or a matrix containing the coefficients for the Finite Element basis expansion. The number of rows
(or the vector's length) corresponds to the number of basis in |
FEMbasis |
A |
An FEM object. This contains a list with components coeff and FEMbasis.
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(horseshoe2D$boundary_nodes, horseshoe2D$locations), segments = horseshoe2D$boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot it plot(FEMfunction)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(horseshoe2D$boundary_nodes, horseshoe2D$locations), segments = horseshoe2D$boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot it plot(FEMfunction)
This function defines a FEM.time object.
FEM.time(coeff,time_mesh,FEMbasis,FLAG_PARABOLIC=FALSE)FEM.time(coeff,time_mesh,FEMbasis,FLAG_PARABOLIC=FALSE)
coeff |
A vector or a matrix containing the coefficients for the spatio-temporal basis expansion. The number of rows
(or the vector's length) corresponds to the number of basis in |
time_mesh |
A vector containing the b-splines knots for separable smoothing and the nodes for finite differences for parabolic smoothing |
FEMbasis |
A |
FLAG_PARABOLIC |
Boolean. If |
A FEM.time object. This contains a list with components coeff, mesh_time, FEMbasis and FLAG_PARABOLIC.
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(horseshoe2D$boundary_nodes, horseshoe2D$locations), segments = horseshoe2D$boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5) time_mesh = seq(0,1,length.out = 5) ## Create the FEM object FEMfunction = FEM.time(coeff, time_mesh, FEMbasis, FLAG_PARABOLIC = TRUE) ## Plot it at desired time plot(FEMfunction,0.7)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(horseshoe2D$boundary_nodes, horseshoe2D$locations), segments = horseshoe2D$boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5) time_mesh = seq(0,1,length.out = 5) ## Create the FEM object FEMfunction = FEM.time(coeff, time_mesh, FEMbasis, FLAG_PARABOLIC = TRUE) ## Plot it at desired time plot(FEMfunction,0.7)
This function implements a smooth functional principal component analysis over a planar mesh, a smooth manifold or a volume.
FPCA.FEM(locations = NULL, datamatrix, FEMbasis, lambda, nPC = 1, validation = NULL, NFolds = 5,GCVmethod = "Stochastic", nrealizations = 100, search = "tree", bary.locations = NULL)FPCA.FEM(locations = NULL, datamatrix, FEMbasis, lambda, nPC = 1, validation = NULL, NFolds = 5,GCVmethod = "Stochastic", nrealizations = 100, search = "tree", bary.locations = NULL)
locations |
A #observations-by-2 matrix in the 2D case and #observations-by-3 matrix in the 2.5D and 3D case, where
each row specifies the spatial coordinates |
datamatrix |
A matrix of dimensions #samples-by-#locations with the observed data values over the domain
for each sample. The datamatrix needs to have zero mean.
If the |
FEMbasis |
A |
lambda |
A scalar or vector of smoothing parameters. |
nPC |
An integer specifying the number of Principal Components to compute. |
validation |
A string specifying the type of validation to perform. If |
NFolds |
This parameter is used only in case |
GCVmethod |
This parameter is considered only when |
nrealizations |
The number of realizations to be used in the stochastic algorithm for the estimation of GCV. |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). |
bary.locations |
A list with three vectors:
|
A list with the following variables:
loadings.FEMA FEM object that represents the L^2-normalized functional loadings for each
Principal Component computed.
scoresA #samples-by-#PrincipalComponents matrix that represents the unnormalized scores or PC vectors.
lambdaA vector of length #PrincipalComponents with the values of the smoothing parameter lambda
chosen for that Principal Component.
variance_explainedA vector of length #PrincipalComponents where each value represent the variance explained by that component.
cumsum_percentageA vector of length #PrincipalComponents containing the cumulative percentage of the variance explained by the first components.
bary.locationsA barycenter information of the given locations if the locations are not mesh nodes.
Lila, E., Aston, J.A.D., Sangalli, L.M., 2016a. Smooth Principal Component Analysis over two-dimensional manifolds with an application to neuroimaging. Ann. Appl. Stat., 10(4), pp. 1854-1879.
library(fdaPDE) ## Load the hub data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) ## Create the Finite Element basis FEMbasis = create.FEM.basis(mesh) ## Create a datamatrix datamatrix = NULL for(ii in 1:50){ a1 = rnorm(1, mean = 1, sd = 1) a2 = rnorm(1, mean = 1, sd = 1) a3 = rnorm(1, mean = 1, sd = 1) func_evaluation = numeric(nrow(mesh$nodes)) for (i in 0:(nrow(mesh$nodes)-1)){ func_evaluation[i+1] = a1* sin(2*pi*mesh$nodes[i+1,1]) + a2* sin(2*pi*mesh$nodes[i+1,2]) + a3*sin(2*pi*mesh$nodes[i+1,3]) + 1 } data = func_evaluation + rnorm(nrow(mesh$nodes), mean = 0, sd = 0.5) datamatrix = rbind(datamatrix, data) } ## Compute the mean of the datamatrix and subtract it to the data data_bar = colMeans(datamatrix) data_demean = matrix(rep(data_bar,50), nrow=50, byrow=TRUE) datamatrix_demeaned = datamatrix - data_demean ## Set the smoothing parameter lambda lambda = 0.00375 ## Estimate the first 2 Principal Components FPCA_solution = FPCA.FEM(datamatrix = datamatrix_demeaned, FEMbasis = FEMbasis, lambda = lambda, nPC = 2) ## Plot the functional loadings of the estimated Principal Components plot(FPCA_solution$loadings.FEM)library(fdaPDE) ## Load the hub data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) ## Create the Finite Element basis FEMbasis = create.FEM.basis(mesh) ## Create a datamatrix datamatrix = NULL for(ii in 1:50){ a1 = rnorm(1, mean = 1, sd = 1) a2 = rnorm(1, mean = 1, sd = 1) a3 = rnorm(1, mean = 1, sd = 1) func_evaluation = numeric(nrow(mesh$nodes)) for (i in 0:(nrow(mesh$nodes)-1)){ func_evaluation[i+1] = a1* sin(2*pi*mesh$nodes[i+1,1]) + a2* sin(2*pi*mesh$nodes[i+1,2]) + a3*sin(2*pi*mesh$nodes[i+1,3]) + 1 } data = func_evaluation + rnorm(nrow(mesh$nodes), mean = 0, sd = 0.5) datamatrix = rbind(datamatrix, data) } ## Compute the mean of the datamatrix and subtract it to the data data_bar = colMeans(datamatrix) data_demean = matrix(rep(data_bar,50), nrow=50, byrow=TRUE) datamatrix_demeaned = datamatrix - data_demean ## Set the smoothing parameter lambda lambda = 0.00375 ## Estimate the first 2 Principal Components FPCA_solution = FPCA.FEM(datamatrix = datamatrix_demeaned, FEMbasis = FEMbasis, lambda = lambda, nPC = 2) ## Plot the functional loadings of the estimated Principal Components plot(FPCA_solution$loadings.FEM)
Implements a finite area test function based on one proposed by Tim Ramsay (2002) proposed by Simon Wood (2008).
fs.test(x, y, r0 = 0.1, r = 0.5, l = 3, b = 1, exclude = FALSE)fs.test(x, y, r0 = 0.1, r = 0.5, l = 3, b = 1, exclude = FALSE)
x, y
|
Points at which to evaluate the test function. |
r0 |
The test domain is a sort of bent sausage. This is the radius of the inner bend. |
r |
The radius of the curve at the centre of the sausage. |
l |
The length of an arm of the sausage. |
b |
The rate at which the function increases per unit increase in distance along the centre line of the sausage. |
exclude |
Should exterior points be set to NA? |
Returns function evaluations, or NAs for points outside the horseshoe domain.
Ramsay, T. 2002. Spline smoothing over difficult regions. J.R.Statist. Soc. B 64(2):307-319
Wood, S. N., Bravington, M. V., & Hedley, S. L. (2008). Soap film smoothing. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(5), 931-955.
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2], exclude = FALSE) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot it plot(FEMfunction)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2], exclude = FALSE) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot it plot(FEMfunction)
Implements a finite area test function based on one proposed by Tim Ramsay (2002) and by Simon Wood (2008) in 3D.
fs.test.3D(x, y, z, r0 = 0.25, r = 1.25, l = 5, b = 1, exclude = FALSE)fs.test.3D(x, y, z, r0 = 0.25, r = 1.25, l = 5, b = 1, exclude = FALSE)
x, y, z
|
Points at which to evaluate the test function. |
r0 |
The test domain is a sort of bent sausage. This is the radius of the inner bend. |
r |
The radius of the curve at the centre of the sausage. |
l |
The length of an arm of the sausage. |
b |
The rate at which the function increases per unit increase in distance along the centre line of the sausage. |
exclude |
Should exterior points be set to NA? |
Returns function evaluations, or NAs for points outside the horseshoe domain.
library(fdaPDE) data(horseshoe2.5D) mesh = horseshoe2.5D FEMbasis=create.FEM.basis(mesh) # Evaluation at nodes sol_exact=fs.test.3D(mesh$nodes[,1],mesh$nodes[,3],mesh$nodes[,2]) plot(FEM(sol_exact, FEMbasis))library(fdaPDE) data(horseshoe2.5D) mesh = horseshoe2.5D FEMbasis=create.FEM.basis(mesh) # Evaluation at nodes sol_exact=fs.test.3D(mesh$nodes[,1],mesh$nodes[,3],mesh$nodes[,2]) plot(FEM(sol_exact, FEMbasis))
A mesh2.5D object with nodes and connectivity matrix of a triangular mesh of the horseshoe 2.5D domain.
The boundary and interior nodes and connectivity matrix of a triangular mesh of the horseshoe domain.
This dataset can be used to create a mesh.2D object with the function create.mesh.2D.
The nodes and connectivity matrix of a triangular mesh of a manifold representing a hub geometry.
This dataset can be used to create a MESH.2.5D object with the function create.MESH.2.5D.
Image plot of a FEM object, generated by the function FEM or returned by
smooth.FEM and FPCA.FEM. Only FEM objects defined over a 2D mesh can be plotted with this method.
## S3 method for class 'FEM' image(x, num_refinements, ...)## S3 method for class 'FEM' image(x, num_refinements, ...)
x |
A 2D-mesh |
num_refinements |
A natural number specifying how many bisections should by applied to each triangular element for plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied. |
... |
Arguments representing graphical options to be passed to plot3d. |
No return value
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot the FEM function image(FEMfunction)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot the FEM function image(FEMfunction)
Image plot of a FEM.time object, generated by the function FEM.time or returned by
smooth.FEM.time. Only FEM objects defined over a 2D mesh can be plotted with this method.
## S3 method for class 'FEM.time' image(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)## S3 method for class 'FEM.time' image(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)
x |
A 2D-mesh |
t |
time at which do the plot |
lambdaS |
index of the space penalization parameter to use for the plot, useful when |
lambdaT |
index of the time penalization parameter to use for the plot, useful when |
num_refinements |
A natural number specifying how many bisections should by applied to each triangular element for plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied. |
... |
Arguments representing graphical options to be passed to plot3d. |
No return value
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 time = 1:5 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time ## Create the FEM.time object FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5,FEMbasis=FEMbasis,FLAG_PARABOLIC=TRUE) ## Plot the FEM function t = c(1.2,1.5,3.6,2.4,4.5) image(FEM_time_function,t)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 time = 1:5 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time ## Create the FEM.time object FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5,FEMbasis=FEMbasis,FLAG_PARABOLIC=TRUE) ## Plot the FEM function t = c(1.2,1.5,3.6,2.4,4.5) image(FEM_time_function,t)
A class that contains all possible information for inference over linear parameters and/or nonparametric field in spatial regression with
differential regularization problem. This object can be used as parameter in smoothing function of the fdaPDE library smooth.FEM.
At least one between test and interval must be nonzero. n_cov, coeff and beta0, if provided, need to be coherent.
dim and locations, if provided, need to be coherent.
The usage of inferenceDataObjectBuilder is recommended for the construction of an object of this class.
testA vector of integers taking value 0, 1 or 2; if 0 no test is performed, if 1 one-at-the-time tests are performed, if 2 a simultaneous test is performed.
intervalA vector of integers taking value 0, 1, 2 or 3; if 0 no confidence interval is computed, if 1 one-at-the-time confidence intervals are computed, if 2 simultaneous confidence intervals are computed, if 3 Bonferroni confidence intervals are computed.
typeA vector of integers taking value 1, 2, 3, 4 or 5 corresponding to Wald, Speckman, Eigen-Sign-Flip, Enhanced-Eigen-Sign-Flip or Sign-Flip inferential approach.
componentA vector of integers taking value 1, 2 or 3, indicating whether the inferential analysis should be carried out respectively for the parametric, nonparametric or both the components.
exactAn integer taking value 1 or 2. If 1 an exact computation of the test statistics will be performed, whereas if 2 an approximated computation will be carried out (not implemented in this version).
dimDimension of the problem, it is equal to 2 in the 1.5D and 2D cases and equal to 3 in the 2.5D and 3D cases.
n_covNumber of covariates taken into account in the linear part of the regression problem.
locationsA matrix of numeric coefficients with columns of dimension dim. When nonparametric inference is requested it represents the set of spatial locations for which the inferential analysis should be performed.
The default values is a one-dimensional matrix of value 1 indicating that all the observed location points should be considered.
In the sign-flip and eigen-sign-flip implementations only observed points are allowed.
locations_indicesA vector of indices indicating which spatial points have to be considered among the observed ones for nonparametric inference. If a vector of location indices is provided then the slot 'location' is discarded.
locations_are_nodesAn integer taking value 1 or 2; in the first case it indicates that the selected locations to perform inference on f are all coinciding with the nodes; otherwise it takes value 2;
coeffA matrix of numeric coefficients with columns of dimension n_cov and each row represents a linear combination of the linear parameters to be tested and/or to be
estimated via confidence interval.
beta0Vector of null hypothesis values for the linear parameters of the model. Used only if test is not 0 and component is not 2.
f0Function representing the expression of the nonparametric component f under the null hypothesis. Used only if component is not 1.
f0_evalVector of f0 evaluations at the chosen test locations. It will be eventually set later in checkInferenceParameters, if nonparametric inference is required.
f_varAn integer taking value 1 or 2. If 1 local variance estimates for the nonlinear part of the model will be computed, whereas if 2 they will not.
quantileVector of quantiles needed for confidence intervals, used only if interval is not 0.
alpha1 minus confidence level vector of sign-flipping approaches confidence intervals. Used only if interval is not 0.
n_flipAn integer representing the number of sign-flips in the case of sign-flipping approaches.
tol_fspaiA real number greater than 0 specifying the tolerance for FSPAI algorithm, in case of non-exact inference (not implemented in this version).
definitionAn integer taking value 0 or 1. If set to 1, the class will be considered as created by the function inferenceDataObjectBuilder,
leading to avoid some of the checks that are performed on inference data within smoothing functions.
A function that build an inferenceDataObject. In the process of construction many checks over the input parameters are carried out so that the output is a well defined object,
that can be used as parameter in smooth.FEM or smooth.FEM.time functions. Notice that this constructor ensures well-posedness of the object, but a further check on consistency with the smoothing functions parameters will be carried out.
inferenceDataObjectBuilder(test = NULL, interval = NULL, type = 'w', component = 'parametric', dim = NULL, n_cov = NULL, locations = NULL, locations_indices = NULL, locations_by_nodes = FALSE, coeff = NULL, beta0 = NULL, f0 = NULL, f_var = FALSE, level = 0.95, n_flip = 1000)inferenceDataObjectBuilder(test = NULL, interval = NULL, type = 'w', component = 'parametric', dim = NULL, n_cov = NULL, locations = NULL, locations_indices = NULL, locations_by_nodes = FALSE, coeff = NULL, beta0 = NULL, f0 = NULL, f_var = FALSE, level = 0.95, n_flip = 1000)
test |
A string defining the type of test to be performed. Multiple tests can be required. In this case the length of the list needs to be coherent with the ones of
|
interval |
A string defining the type of confidence intervals to be computed. Multiple intervals can be required. In this case the length of the list needs to be coherent with the ones of
|
type |
A list of strings defining the type of implementation for the inferential analysis. The possible values are:
|
component |
A list of strings defining on which model component inference has to be performed. It can take values 'parametric' (default), 'nonparametric' or 'both'. |
dim |
Dimension of the problem, defaulted to NULL. It can take value 2 or 3 corresponding to 1.5D/2D or 2.5D/3D problems (Must be set by the user) |
n_cov |
Number of the covariates, defaulted to NULL. (Must be set by the user) |
locations |
A matrix of the locations of interest when testing the nonparametric component f, defaulted to NULL |
locations_indices |
A vector of indices indicating the locations to be considered among the observed ones for nonparametric inference, defaulted to NULL. If a vector of indices is provided, then the slot 'locations' is discarded. |
locations_by_nodes |
A logical used to indicate whether the selected locations to perform inference on f are all coinciding with the nodes; |
coeff |
A matrix, with |
beta0 |
Vector of real numbers (default NULL). It is used only if the |
f0 |
A function object representing the expression of the nonparametric component f under the null hypothesis. Taken into account if |
f_var |
A logical used to decide whether to estimate the local variance of the nonlinear part of the model. The possible values are: FALSE (default) and TRUE. |
level |
A vector containing the level of significance used to compute quantiles for confidence intervals, defaulted to 0.95. It is taken into account only if |
n_flip |
Number of flips performed in sign-flipping approaches, defaulted to 1000. |
The output is a well defined inferenceDataObject, that can be used as input parameter in the smooth.FEM function.
obj<-inferenceDataObjectBuilder(test = 'oat', dim = 2, beta0 = rep(1,4), n_cov = 4); obj2<-inferenceDataObjectBuilder(test = 'sim', dim = 3, component = 'nonparametric', n_cov = 3);obj<-inferenceDataObjectBuilder(test = 'oat', dim = 2, beta0 = rep(1,4), n_cov = 4); obj2<-inferenceDataObjectBuilder(test = 'sim', dim = 3, component = 'nonparametric', n_cov = 3);
A class that contains all possible information for inference over linear parameters and/or nonparametric field in spatio-temporal regression with
differential regularization problem. This object can be used as parameter in smoothing function of the fdaPDE library smooth.FEM.time.
At least one between test and interval must be nonzero. n_cov, coeff and beta0, if provided, need to be coherent.
dim and locations, if provided, need to be coherent.
The usage of inferenceDataObjectTimeBuilder is recommended for the construction of an object of this class.
testA vector of integers taking value 0, 1 or 2; if 0 no test is performed, if 1 one-at-the-time tests are performed, if 2 a simultaneous test is performed.
intervalA vector of integers taking value 0, 1, 2 or 3; if 0 no confidence interval is computed, if 1 one-at-the-time confidence intervals are computed, if 2 simultaneous confidence intervals are computed, if 3 Bonferroni confidence intervals are computed.
typeA vector of integers taking value 1, 2, 3 or 4 corresponding to Wald, Speckman, Eigen-Sign-Flip, Enhanced-Eigen-Sign-Flip inferential approach.
componentA vector of integers taking value 1, 2 or 3, indicating whether the inferential analysis should be carried out respectively for the parametric, nonparametric or both the components.
exactAn integer taking value 1 or 2. If 1 an exact computation of the test statistics will be performed, whereas if 2 an approximated computation will be carried out (not implemented in this version).
dimDimension of the problem, it is equal to 2 in the 1.5D and 2D cases and equal to 3 in the 2.5D and 3D cases.
n_covNumber of covariates taken into account in the linear part of the regression problem.
locationsA matrix of numeric coefficients with columns of dimension dim. When nonparametric inference is requested it represents the set of spatial locations for which the inferential analysis should be performed.
The default values is a one-dimensional matrix of value 1 indicating that all the observed location points should be considered.
In the sign-flip and eigen-sign-flip implementations only observed points are allowed.
locations_indicesA vector of indices indicating which spatial points have to be considered among the observed ones for nonparametric inference. If a vector of location indices is provided then the slot 'location' is discarded.
locations_are_nodesAn integer taking value 1 or 2; in the first case it indicates that the selected locations to perform inference on f are all coinciding with the nodes; otherwise it takes value 2;
time_locationsA vector of numeric coefficients containing the set of times of interest for inference on the nonparametric component. This parameter can be NULL.
In this case the temporal locations are assumed to coincide with the time_locations provided to the smoothing functions. Used only if component is not 1.
coeffA matrix of numeric coefficients with columns of dimension n_cov and each row represents a linear combination of the linear parameters to be tested and/or to be
estimated via confidence interval.
beta0Vector of null hypothesis values for the linear parameters of the model. Used only if test is not 0 and component is not 2.
f0Function representing the expression of the nonparametric component f under the null hypothesis. Used only if component is not 1.
f0_evalMatrix of f0 evaluations at the chosen space and time locations. It will be eventually set later in checkInferenceParametersTime, if nonparametric inference is required.
f_varAn integer taking value 1 or 2. If 1 local variance estimates for the nonlinear part of the model will be computed, whereas if 2 they will not.
quantileVector of quantiles needed for confidence intervals, used only if interval is not 0.
alpha1 minus confidence level vector of sign-flipping approaches confidence intervals. Used only if interval is not 0.
n_flipAn integer representing the number of sign-flips in the case of sign-flipping approaches.
tol_fspaiA real number greater than 0 specifying the tolerance for FSPAI algorithm, in case of non-exact inference (not implemented in this version).
definitionAn integer taking value 0 or 1. If set to 1, the class will be considered as created by the function inferenceDataObjectTimeBuilder,
leading to avoid some of the checks that are performed on inference data within smoothing functions.
A function that build an inferenceDataObjectTime. In the process of construction many checks over the input parameters are carried out so that the output is a well defined object,
that can be used as parameter in smooth.FEM or smooth.FEM.time functions.
Notice that this constructor ensures well-posedness of the object, but a further check on consistency with the smoothing functions parameters will be carried out.
inferenceDataObjectTimeBuilder(test = NULL, interval = NULL, type = 'w', component = 'parametric', dim = NULL, n_cov = NULL, locations = NULL, locations_indices = NULL, locations_by_nodes = F, time_locations = NULL, coeff = NULL, beta0 = NULL, f0 = NULL, f_var = F, level = 0.95, n_flip = 1000)inferenceDataObjectTimeBuilder(test = NULL, interval = NULL, type = 'w', component = 'parametric', dim = NULL, n_cov = NULL, locations = NULL, locations_indices = NULL, locations_by_nodes = F, time_locations = NULL, coeff = NULL, beta0 = NULL, f0 = NULL, f_var = F, level = 0.95, n_flip = 1000)
test |
A string defining the type of test to be performed. Multiple tests can be required. In this case the length of the list needs to be coherent with the ones of
|
interval |
A string defining the type of confidence intervals to be computed. Multiple intervals can be required. In this case the length of the list needs to be coherent with the ones of
|
type |
A list of strings defining the type of implementation for the inferential analysis. The possible values are:
|
component |
A list of strings defining on which model component inference has to be performed. It can take values 'parametric' (default), 'nonparametric' or 'both'. |
dim |
Dimension of the problem, defaulted to NULL. It can take value 2 or 3 corresponding to 1.5D/2D or 2.5D/3D problems (Must be set by the user) |
n_cov |
Number of the covariates, defaulted to NULL. (Must be set by the user) |
locations |
A matrix of the locations of interest when testing the nonparametric component f, defaulted to NULL. |
locations_indices |
A vector of indices indicating the locations to be considered among the observed ones for nonparametric inference, defaulted to NULL. If a vector of indices is provided, then the slot 'locations' is discarded. |
locations_by_nodes |
A logical used to indicate whether the selected locations to perform inference on f are all coinciding with the nodes. |
time_locations |
A vector of times of interest when testing the nonparametric component f, defaulted to NULL. If |
coeff |
A matrix, with |
beta0 |
Vector of real numbers (default NULL). It is used only if the |
f0 |
A function object representing the expression of the nonparametric component f under the null hypothesis. Taken into account if |
f_var |
A logical used to decide whether to estimate the local variance of the nonlinear part of the model. The possible values are: FALSE (default) and TRUE. |
level |
A vector containing the level of significance used to compute quantiles for confidence intervals, defaulted to 0.95. It is taken into account only if |
n_flip |
Number of flips performed in sign-flipping approaches, defaulted to 1000. |
The output is a well defined inferenceDataObjectTime, that can be used as input parameter in the smooth.FEM function.
obj<-inferenceDataObjectTimeBuilder(test = 'oat', dim = 2, beta0 = rep(1,4), n_cov = 4); obj2<-inferenceDataObjectTimeBuilder(test = 'sim', dim = 3, component = 'nonparametric', n_cov = 3);obj<-inferenceDataObjectTimeBuilder(test = 'oat', dim = 2, beta0 = rep(1,4), n_cov = 4); obj2<-inferenceDataObjectTimeBuilder(test = 'sim', dim = 3, component = 'nonparametric', n_cov = 3);
FEM objectThree-dimensional plot of a FEM object, generated by FEM or returned by
smooth.FEM or FPCA.FEM.
If the mesh of the FEMbasis component is of class mesh.2D both the 3rd axis and the color represent
the value of the coefficients for the Finite Element basis expansion (coeff component of the FEM object).
If the mesh is of class mesh.3D, the color of each triangle or tetrahedron represent the mean value of
the coefficients for the Finite Element basis expansion (coeff).
## S3 method for class 'FEM' plot(x, colormap = "heat.colors", num_refinements = NULL, ...)## S3 method for class 'FEM' plot(x, colormap = "heat.colors", num_refinements = NULL, ...)
x |
A |
colormap |
A colormap exploited in the plot. The default value is the heat colormap. |
num_refinements |
A natural number specifying how many bisections should be applied to each triangular element for plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied. This parameter can be specified only when a FEM object defined over a 2D mesh is plotted. |
... |
Arguments representing graphical options to be passed to plot3d. |
No return value
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot the FEM function plot(FEMfunction)library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2]) ## Create the FEM object FEMfunction = FEM(coeff, FEMbasis) ## Plot the FEM function plot(FEMfunction)
FEM.time object at a given timePlot of a FEM.time object, generated by FEM.time or returned by
smooth.FEM.time. time_locations and locations must not be both provided.
If time_locations is provided, the spatial field is plotted for the provided temporal instnts.
If locations is provided, the temporal evolution in the provided space locations is plotted.
If both time_locations and locations are NULL a default plot is provided.
If the mesh of the FEMbasis component is of class mesh.2D both the 3rd axis and the color represent
the value of the coefficients for the Finite Element basis expansion (coeff component of the FEM.time object).
If the mesh is of class mesh.3D, the color of each triangle or tetrahedron represent the mean value of
the coefficients for the Finite Element basis expansion (coeff).
## S3 method for class 'FEM.time' plot(x, time_locations = NULL, locations = NULL, lambdaS = NULL, lambdaT = NULL, num_refinements = NULL, Nt = 100, add = FALSE, main = NULL, col = "red", ...)## S3 method for class 'FEM.time' plot(x, time_locations = NULL, locations = NULL, lambdaS = NULL, lambdaT = NULL, num_refinements = NULL, Nt = 100, add = FALSE, main = NULL, col = "red", ...)
x |
A |
time_locations |
A vector with the instants in which plot the spatial field |
locations |
A 2-column matrix when |
lambdaS |
Index of the space penalization parameter to use for the plot, useful when |
lambdaT |
Index of the time penalization parameter to use for the plot, useful when |
num_refinements |
A natural number specifying how many bisections should be applied to each triangular element for plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied. This parameter can be specified only when a FEM object defined over a 2D mesh is plotted. |
Nt |
The number of instants to plot when |
add |
Boolean, used only when locations is not NULL, if TURE it performs the graphic over the last plot |
main |
The title of the plot when |
col |
The color of the plot when |
... |
Arguments representing graphical options to be passed to plot3d. |
No return value
library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 time = 1:5 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time ## Create the FEM.time object FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5, FEMbasis=FEMbasis, FLAG_PARABOLIC=TRUE) ## Plot the FEM function plot(FEM_time_function) ## plot spatial field in some instants t = c(1.2,1.5,3.6,2.4,4.5) plot(FEM_time_function, t) ## plot time evolution in some locations plot(FEM_time_function, locations = locations[1:10,])library(fdaPDE) ## Upload the horseshoe2D data data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations ## Create the 2D mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Create the FEM basis FEMbasis = create.FEM.basis(mesh) ## Compute the coeff vector evaluating the desired function at the mesh nodes ## In this case we consider the fs.test() function introduced by Wood et al. 2008 time = 1:5 coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time ## Create the FEM.time object FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5, FEMbasis=FEMbasis, FLAG_PARABOLIC=TRUE) ## Plot the FEM function plot(FEM_time_function) ## plot spatial field in some instants t = c(1.2,1.5,3.6,2.4,4.5) plot(FEM_time_function, t) ## plot time evolution in some locations plot(FEM_time_function, locations = locations[1:10,])
Plot a mesh.1.5D object, generated by create.mesh.1.5D or refine.mesh.1.5D.
## S3 method for class 'mesh.1.5D' plot(x, ...)## S3 method for class 'mesh.1.5D' plot(x, ...)
x |
A |
... |
Arguments representing graphical options to be passed to par. |
No return value
Plot the triangulation of a mesh.2.5D object, generated by create.mesh.2.5D
## S3 method for class 'mesh.2.5D' plot(x, ...)## S3 method for class 'mesh.2.5D' plot(x, ...)
x |
A |
... |
Arguments representing graphical options to be passed to par. |
No return value
library(fdaPDE) ## Upload the hub2.5D the data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create mesh mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) plot(mesh)library(fdaPDE) ## Upload the hub2.5D the data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create mesh mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) plot(mesh)
Plot a mesh.2D object, generated by create.mesh.2D or refine.mesh.2D.
## S3 method for class 'mesh.2D' plot(x, ...)## S3 method for class 'mesh.2D' plot(x, ...)
x |
A |
... |
Arguments representing graphical options to be passed to par. |
No return value
library(fdaPDE) ## Upload the quasicirle2D data data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data ## Create mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Plot the mesh plot(mesh)library(fdaPDE) ## Upload the quasicirle2D data data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data ## Create mesh mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) ## Plot the mesh plot(mesh)
Plot a mesh.3D object, generated by create.mesh.3D.
## S3 method for class 'mesh.3D' plot(x, ...)## S3 method for class 'mesh.3D' plot(x, ...)
x |
A |
... |
Arguments representing graphical options to be passed to par. |
No return value
library(fdaPDE) ##Load the matrix nodes and tetrahedrons data(sphere3Ddata) nodes = sphere3Ddata$nodes tetrahedrons = sphere3Ddata$tetrahedrons ##Create the triangulated mesh from the connectivity matrix and nodes locations mesh = create.mesh.3D(nodes,tetrahedrons) ##Plot the triangulation of the object plot(mesh)library(fdaPDE) ##Load the matrix nodes and tetrahedrons data(sphere3Ddata) nodes = sphere3Ddata$nodes tetrahedrons = sphere3Ddata$tetrahedrons ##Create the triangulated mesh from the connectivity matrix and nodes locations mesh = create.mesh.3D(nodes,tetrahedrons) ##Plot the triangulation of the object plot(mesh)
This function projects any 2D points onto 1.5D linear network mesh.
projection.points.1.5D(mesh, locations)projection.points.1.5D(mesh, locations)
mesh |
A mesh.1.5D object representing the graph mesh, created by create.mesh.1.5D. |
locations |
2D points to be projected onto 1.5D mesh. |
2D points projected onto 1.5D linear network mesh.
library(fdaPDE) ##Create Mesh nodes=matrix(c(0.25,0.25,0.5,0.25,0.75,0.5,0.75,0.), nrow = 4, byrow=TRUE) edges=matrix(c(1,2,2,3,2,4),nrow = 3,byrow = TRUE) mesh_ = create.mesh.1.5D(nodes,edges,order=1) ## Create 2D points to be projected locations=matrix(nrow=5,ncol=2) locations[,1] = runif(5,min=0.25,max=0.75) locations[,2] = runif(5,min=0.25,max=0.5) ## Project the points on the mesh loc = projection.points.1.5D(mesh_, locations)library(fdaPDE) ##Create Mesh nodes=matrix(c(0.25,0.25,0.5,0.25,0.75,0.5,0.75,0.), nrow = 4, byrow=TRUE) edges=matrix(c(1,2,2,3,2,4),nrow = 3,byrow = TRUE) mesh_ = create.mesh.1.5D(nodes,edges,order=1) ## Create 2D points to be projected locations=matrix(nrow=5,ncol=2) locations[,1] = runif(5,min=0.25,max=0.75) locations[,2] = runif(5,min=0.25,max=0.5) ## Project the points on the mesh loc = projection.points.1.5D(mesh_, locations)
This function projects any 3D points onto 2.5D triangular mesh.
projection.points.2.5D(mesh, locations)projection.points.2.5D(mesh, locations)
mesh |
A mesh.2.5D object representing the triangular mesh, created by create.mesh.2.5D. |
locations |
3D points to be projected onto 2.5D triangular mesh. |
3D points projected onto 2.5D triangluar mesh.
library(fdaPDE) ## Upload the hub2.5D the data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create mesh mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) ## Create 3D points to be projected x <- cos(seq(0,2*pi, length.out = 9)) y <- sin(seq(0,2*pi, length.out = 9)) z <- rep(0.5, 9) locations = cbind(x,y,z) ## Project the points on the mesh loc = projection.points.2.5D(mesh, locations)library(fdaPDE) ## Upload the hub2.5D the data data(hub2.5D) hub2.5D.nodes = hub2.5D$hub2.5D.nodes hub2.5D.triangles = hub2.5D$hub2.5D.triangles ## Create mesh mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles) ## Create 3D points to be projected x <- cos(seq(0,2*pi, length.out = 9)) y <- sin(seq(0,2*pi, length.out = 9)) z <- rep(0.5, 9) locations = cbind(x,y,z) ## Project the points on the mesh loc = projection.points.2.5D(mesh, locations)
The boundary and interior nodes and connectivity matrix of a triangular mesh of a quasicircular domain, together
with a non-stationary field observed over the nodes of the mesh.
This dataset can be used to create a mesh.2D object with the function create.mesh.2D and to test
the smooth.FEM function.
The mesh of a quasicircular domain, together with a non-stationary field observed over seven circular subdomains and the incindence matrix defining the subdomains used by Azzimonti et. al 2015. This dataset can be used to test the smooth.FEM function for areal data.
Azzimonti, L., Sangalli, L. M., Secchi, P., Domanin, M., & Nobile, F. (2015). Blood flow velocity field estimation via spatial regression with PDE penalization. Journal of the American Statistical Association, 110(511), 1057-1071.
mesh.1.5D object by splitting each edge of a given mesh into two subedges.Create a mesh.1.5D object by splitting each edge of a given mesh into two subedges.
refine.by.splitting.mesh.1.5D(mesh = NULL)refine.by.splitting.mesh.1.5D(mesh = NULL)
mesh |
a |
An object of class mesh.1.5D with splitted edges
mesh.2.5D object by splitting each triangle of a given mesh into four subtriangles.Create a mesh.2.5D object by splitting each triangle of a given mesh into four subtriangles.
refine.by.splitting.mesh.2.5D(mesh = NULL)refine.by.splitting.mesh.2.5D(mesh = NULL)
mesh |
a |
An object of class mesh.2.5D with splitted triangles
mesh.2D object by splitting each triangle of a given mesh into four subtriangles.Create a mesh.2D object by splitting each triangle of a given mesh into four subtriangles.
refine.by.splitting.mesh.2D(mesh = NULL)refine.by.splitting.mesh.2D(mesh = NULL)
mesh |
a |
An object of class mesh.2D with splitted triangles
mesh.3D object by splitting each tetrahedron of a given mesh into eight subtetrahedrons.Create a mesh.3D object by splitting each tetrahedron of a given mesh into eight subtetrahedrons.
refine.by.splitting.mesh.3D(mesh = NULL)refine.by.splitting.mesh.3D(mesh = NULL)
mesh |
a |
An object of class mesh.3D with splitted tetrahedrons
Refine 1.5D mesh
refine.mesh.1.5D(mesh, delta)refine.mesh.1.5D(mesh, delta)
mesh |
a |
delta |
the maximum allowed length |
An object of class mesh.1.5D with refined edges
This function refines a Constrained Delaunay triangulation into a Conforming Delaunay triangulation. This is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used to
refine a mesh previously created with create.mesh.2D. The algorithm can add Steiner points (points through which the segments are splitted)
in order to meet the imposed refinement conditions.
refine.mesh.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity)refine.mesh.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity)
mesh |
A mesh.2D object representing the triangular mesh, created by create.mesh.2D. |
minimum_angle |
A scalar specifying a minimun value for the triangles angles. |
maximum_area |
A scalar specifying a maximum value for the triangles areas. |
delaunay |
A boolean parameter indicating whether or not the output mesh should satisfy the Delaunay condition. |
verbosity |
This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process. |
A mesh.2D object representing the refined triangular mesh, with the following output:
nodesA #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
nodesmarkersA vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.
nodesattributesnodesattributes A matrix with #nodes rows containing nodes' attributes. These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value of a Dirichlet boundary condition at boundary nodes added during the triangulation process.
trianglesA #triangles-by-3 (when order = 1) or #triangles-by-6 (when order = 2) matrix.
edgesA #edges-by-2 matrix. Each row contains the row's indices of the nodes where the edge starts from and ends to.
edgesmarkersA vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in edge is a boundary edge;
an entry '0' indicates that the corresponding edge is not a boundary edge.
neighborsA #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that one edge of the triangle is a boundary edge.
holesA #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes in the triangulation, when the domain has holes.
orderEither '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints). These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.
create.mesh.2D, create.FEM.basis
library(fdaPDE) ## Upload the quasicircle2D data data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data ## Create mesh from boundary: mesh = create.mesh.2D(nodes = boundary_nodes, segments = boundary_segments) plot(mesh) ## Refine the mesh with the maximum area criterion: finemesh = refine.mesh.2D(mesh = mesh, maximum_area = 0.1) plot(finemesh) ## Refine the mesh with the minimum angle criterion: finemesh2 = refine.mesh.2D(mesh = mesh, minimum_angle = 30) plot(finemesh2)library(fdaPDE) ## Upload the quasicircle2D data data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data ## Create mesh from boundary: mesh = create.mesh.2D(nodes = boundary_nodes, segments = boundary_segments) plot(mesh) ## Refine the mesh with the maximum area criterion: finemesh = refine.mesh.2D(mesh = mesh, maximum_area = 0.1) plot(finemesh) ## Refine the mesh with the minimum angle criterion: finemesh2 = refine.mesh.2D(mesh = mesh, minimum_angle = 30) plot(finemesh2)
This function implements a spatial regression model with differential regularization. The regularizing term involves a Partial Differential Equation (PDE). In the simplest case the PDE involves only the Laplacian of the spatial field, that induces an isotropic smoothing. When prior information about the anisotropy or non-stationarity is available the PDE involves a general second order linear differential operator with possibly space-varying coefficients. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries
smooth.FEM(locations = NULL, observations, FEMbasis, covariates = NULL, PDE_parameters = NULL, BC = NULL, incidence_matrix = NULL, areal.data.avg = TRUE, search = "tree", bary.locations = NULL, family = "gaussian", mu0 = NULL, scale.param = NULL, threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15, lambda.selection.criterion = "grid", DOF.evaluation = NULL, lambda.selection.lossfunction = NULL, lambda = NULL, DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0, DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05, inference.data.object=NULL)smooth.FEM(locations = NULL, observations, FEMbasis, covariates = NULL, PDE_parameters = NULL, BC = NULL, incidence_matrix = NULL, areal.data.avg = TRUE, search = "tree", bary.locations = NULL, family = "gaussian", mu0 = NULL, scale.param = NULL, threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15, lambda.selection.criterion = "grid", DOF.evaluation = NULL, lambda.selection.lossfunction = NULL, lambda = NULL, DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0, DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05, inference.data.object=NULL)
locations |
A #observations-by-2 matrix in the 2D case and #observations-by-3 matrix in the 2.5D and 3D case, where
each row specifies the spatial coordinates |
observations |
A vector of length #observations with the observed data values over the domain.
If the |
FEMbasis |
A |
covariates |
A #observations-by-#covariates matrix where each row represents the covariates associated with
the corresponding observed data value in |
PDE_parameters |
A list specifying the parameters of the PDE in the regularizing term. Default is NULL, i.e.
regularization is by means of the Laplacian (stationary, isotropic case).
If the coefficients of the PDE are constant over the domain
If the coefficients of the PDE are space-varying
For 2.5D and 3D, only the Laplacian is available ( |
BC |
A list with two vectors:
|
incidence_matrix |
A #regions-by-#triangles/tetrahedrons matrix where the element (i,j) equals 1 if the j-th
triangle/tetrahedron is in the i-th region and 0 otherwise.
This is needed only for areal data. In case of pointwise data, this parameter is set to |
areal.data.avg |
Boolean. It involves the computation of Areal Data. If |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). |
bary.locations |
A list with three vectors:
|
family |
This parameter specify the distibution within exponential family used for GLM model.
The following distribution are implemented: "binomial", "exponential", "gamma", "poisson", "gaussian", "invgaussian".
The default link function for binomial is |
mu0 |
This parameter is a vector that set the starting point for FPIRLS algorithm. It represent an initial guess of the location parameter.
Default is set to observation for non binary distribution while equal to |
scale.param |
Dispersion parameter of the chosen distribution. This is only required for "gamma", "gaussian", "invgaussian". User may specify the parameter as a positive real number. If the parameter is not supplied, it is estimated from data according to Wilhelm Sangalli 2016. |
threshold.FPIRLS |
This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.FPIRLS.
Default value |
max.steps.FPIRLS |
This parameter is used to limit the maximum number of iteration.
Default value |
lambda.selection.criterion |
This parameter is used to select the optimization method for the smoothing parameter |
DOF.evaluation |
This parameter is used to identify if and how to perform degrees of freedom computation.
The following possibilities are allowed: NULL, 'exact' and 'stochastic'
In the former case no degree of freedom is computed, while the other two methods enable computation.
Stochastic computation of DOFs may be slightly less accurate than its deterministic counterpart, but it is fairly less time consuming. Stochastic evaluation is highly suggested for meshes with more than 5000 nodes.
Default value |
lambda.selection.lossfunction |
This parameter is used to determine if some loss function has to be evaluated.
The following possibilities are allowed: NULL and 'GCV' (generalized cross validation)
If NULL is selected, |
lambda |
a vector of spatial smoothing parameters provided if |
DOF.stochastic.realizations |
This positive integer is considered only when |
DOF.stochastic.seed |
This positive integer is considered only when |
DOF.matrix |
Matrix of degrees of freedom. This parameter can be used if the DOF.matrix corresponding to |
GCV.inflation.factor |
Tuning parameter used for the estimation of GCV. Default value |
lambda.optimization.tolerance |
Tolerance parameter, a double between 0 and 1 that fixes how much precision is required by the optimization method: the smaller the parameter, the higher the accuracy.
Used only if |
inference.data.object |
An |
A list with the following variables:
fit.FEMA FEM object that represents the fitted spatial field.
PDEmisfit.FEMA FEM object that represents the Laplacian of the estimated spatial field.
solutionA list, note that all terms are matrices or row vectors: the jth column represents the vector related to lambda[j] if lambda.selection.criterion="grid" and lambda.selection.lossfunction=NULL.
In all the other cases, only the column related to the best smoothing parameter is returned.
fMatrix, estimate of function f, first half of solution vector.
gMatrix, second half of solution vector.
z_hatMatrix, prediction of the output in the locations.
betaIf covariates is not NULL, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. It is the regression coefficients estimate.
rmseEstimate of the root mean square error in the locations.
estimated_sdEstimate of the standard deviation of the error.
optimizationA detailed list of optimization related data:
lambda_solutionnumerical value of best lambda according to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction=NULL.
lambda_positioninteger, position in lambda_vector of best lambda according to
lambda.selection.lossfunction, -1 if lambda.selection.lossfunction=NULL.
GCVnumeric value of GCV in correspondence of the optimum.
optimization_detailslist containing further information about the optimization method used and the nature of its termination, eventual number of iterations.
dofvector of positive numbers, DOFs for all the lambdas in lambda_vector, empty or invalid if not computed.
lambda_vectorvector of positive numbers: penalizations either passed by the user or found in the iterations of the optimization method.
GCV_vectorvector of positive numbers, GCV values for all the lambdas in lambda_vector
timeDuration of the entire optimization computation.
bary.locationsBarycenter information of the given locations, if the locations are not mesh nodes.
GAM_outputA list of GAM related data:
fn_hatA matrix with number of rows equal to number of locations and number of columns equal to length of lambda. Each column contains the evaluaton of the spatial field in the location points.
J_minimaA vector of the same length of lambda, containing the reached minima for each value of the smoothing parameter.
variance.estA vector which returns the variance estimates for the Generative Additive Models.
inferenceA list set only if a well defined inferenceDataObject is passed as parameter to the function; contains all inference outputs required:
p_valueslist of lists set only if at least one p-value is required; contains the p-values divided by implementation:
waldlist containing all the Wald p-values required, in the same order of the type list in inference.data.object.
If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of coeff matrix in inference.data.object.
speckmanlist containing all the Speckman p-values required, in the same order of the type list in inference.data.object.
If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of coeff matrix in inference.data.object.
eigen_sign_fliplist containing all the Eigen-Sign-Flip p-values required, in the same order of the type list in inference.data.object.
If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of coeff matrix in inference.data.object.
CIlist of lists set only if at least one confidence interval is required; contains the confidence intervals divided by implementation:
waldlist containing all the Wald confidence intervals required, in the same order of the type list in inference.data.object.
Each item is a matrix with 3 columns and p rows, p being the number of rows of coeff matrix in inference.data.object; each row is the CI for the corresponding row of coeff matrix.
speckmanlist containing all the Speckman confidence intervals required, in the same order of the type list in inference.data.object.
Each item is a matrix with 3 columns and p rows, p being the number of rows of coeff matrix in inference.data.object; each row is the CI for the corresponding row of coeff matrix.
Sangalli, L. M., Ramsay, J. O., Ramsay, T. O. (2013). Spatial spline regression models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(4), 681-703.
Azzimonti, L., Sangalli, L. M., Secchi, P., Domanin, M., Nobile, F. (2015). Blood flow velocity field estimation via spatial regression with PDE penalization. Journal of the American Statistical Association, 110(511), 1057-1071.
Matthieu Wilhelm & Laura M. Sangalli (2016). Generalized spatial regression with differential regularization. Journal of Statistical Computation and Simulation, 86:13, 2497-2518.
Federico Ferraccioli, Laura M. Sangalli & Livio Finos (2022). Some first inferential tools for spatial regression with differential regularization. Journal of Multivariate Analysis, 189, 104866.
library(fdaPDE) #### No prior information about anysotropy/non-stationarity (laplacian smoothing) #### data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) FEMbasis = create.FEM.basis(mesh) lambda = 10^-1 # no covariate data = fs.test(mesh$nodes[,1], mesh$nodes[,2]) + rnorm(nrow(mesh$nodes), sd = 0.5) solution = smooth.FEM(observations = data, FEMbasis = FEMbasis, lambda = lambda) plot(solution$fit.FEM) # with covariates covariate = covs.test(mesh$nodes[,1], mesh$nodes[,2]) data = fs.test(mesh$nodes[,1], mesh$nodes[,2]) + 2*covariate + rnorm(nrow(mesh$nodes), sd = 0.5) #Inferential tests and confidence intervals inference.data.object = inferenceDataObjectBuilder(test = 'oat', type = 'w', dim = 2, n_cov = 1) solution = smooth.FEM(observations = data, covariates = covariate, FEMbasis = FEMbasis, lambda = lambda, inference.data.object=inference.data.object) # beta estimate: solution$solution$beta # tests over beta estimates p-values: solution$inference$beta$p_values # confidence intervals for beta estimates: solution$inference$beta$CI # non-parametric estimate: plot(solution$fit.FEM) # Choose lambda with GCV - stochastic grid evaluation: lambda = 10^(-2:0) solution = smooth.FEM(observations = data, covariates = covariate, FEMbasis = FEMbasis, lambda = lambda, DOF.evaluation = 'stochastic', lambda.selection.lossfunction = 'GCV') bestLambda = solution$optimization$lambda_solution # Choose lambda with GCV - Newton finite differences stochastic evaluation -: solution = smooth.FEM(observations = data, covariates = covariate, FEMbasis = FEMbasis, DOF.evaluation = 'stochastic', lambda.selection.lossfunction = 'GCV') bestLambda = solution$optimization$lambda_solution #### Smoothing with prior information about anysotropy/non-stationarity and boundary conditions #### # See Azzimonti et al. for reference to the current example data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) FEMbasis = create.FEM.basis(mesh) lambda = 10^-2 # Set the PDE parameters R = 2.8 K1 = 0.1 K2 = 0.2 beta = 0.5 K_func<-function(points) { output = array(0, c(2, 2, nrow(points))) for (i in 1:nrow(points)) output[,,i]=10*rbind(c(points[i,2]^2+K1*points[i,1]^2+K2*(R^2-points[i,1]^2-points[i,2]^2), (K1-1)*points[i,1]*points[i,2]), c((K1-1)*points[i,1]*points[i,2], points[i,1]^2+K1*points[i,2]^2+K2*(R^2-points[i,1]^2-points[i,2]^2))) output } b_func<-function(points) { output = array(0, c(2, nrow(points))) for (i in 1:nrow(points)) output[,i] = 10*beta*c(points[i,1],points[i,2]) output } c_func<-function(points) { rep(c(0), nrow(points)) } u_func<-function(points) { rep(c(0), nrow(points)) } PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func) # Set the boundary conditions BC = NULL BC$BC_indices = which(mesh$nodesmarkers == 1) # b.c. on the complete boundary BC$BC_values = rep(0,length(BC$BC_indices)) # homogeneus b.c. # Since the data locations are a subset of the mesh nodes for a faster solution use: dataNA = rep(NA, FEMbasis$nbasis) dataNA[mesh$nodesmarkers == 0] = data #grid evaluation solution = smooth.FEM(observations = dataNA, FEMbasis = FEMbasis, lambda = lambda, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM) # Newton's method solution = smooth.FEM(observations = dataNA, FEMbasis = FEMbasis, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM) #### Smoothing with areal data #### # See Azzimonti et al. for reference to the current exemple data(quasicircle2Dareal) incidence_matrix = quasicircle2Dareal$incidence_matrix data = quasicircle2Dareal$data mesh = quasicircle2Dareal$mesh FEMbasis = create.FEM.basis(mesh) lambda = 10^-4 # Set the PDE parameters R = 2.8 K1 = 0.1 K2 = 0.2 beta = 0.5 K_func<-function(points) { output = array(0, c(2, 2, nrow(points))) for (i in 1:nrow(points)) output[,,i]=10*rbind(c(points[i,2]^2+K1*points[i,1]^2+K2*(R^2-points[i,1]^2-points[i,2]^2), (K1-1)*points[i,1]*points[i,2]), c((K1-1)*points[i,1]*points[i,2], points[i,1]^2+K1*points[i,2]^2+K2*(R^2-points[i,1]^2-points[i,2]^2))) output } b_func<-function(points) { output = array(0, c(2, nrow(points))) for (i in 1:nrow(points)) output[,i] = 10*beta*c(points[i,1],points[i,2]) output } c_func<-function(points) { rep(c(0), nrow(points)) } u_func<-function(points) { rep(c(0), nrow(points)) } PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func) # Set the boundary conditions BC = NULL BC$BC_indices = which(mesh$nodesmarkers == 1) # b.c. on the complete boundary BC$BC_values = rep(0,length(BC$BC_indices)) # homogeneus b.c. #grid evaluation solution = smooth.FEM(observations = data, incidence_matrix = incidence_matrix, FEMbasis = FEMbasis, lambda = lambda, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM) #Newton's method solution = smooth.FEM(observations = data, incidence_matrix = incidence_matrix, FEMbasis = FEMbasis, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM)library(fdaPDE) #### No prior information about anysotropy/non-stationarity (laplacian smoothing) #### data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) FEMbasis = create.FEM.basis(mesh) lambda = 10^-1 # no covariate data = fs.test(mesh$nodes[,1], mesh$nodes[,2]) + rnorm(nrow(mesh$nodes), sd = 0.5) solution = smooth.FEM(observations = data, FEMbasis = FEMbasis, lambda = lambda) plot(solution$fit.FEM) # with covariates covariate = covs.test(mesh$nodes[,1], mesh$nodes[,2]) data = fs.test(mesh$nodes[,1], mesh$nodes[,2]) + 2*covariate + rnorm(nrow(mesh$nodes), sd = 0.5) #Inferential tests and confidence intervals inference.data.object = inferenceDataObjectBuilder(test = 'oat', type = 'w', dim = 2, n_cov = 1) solution = smooth.FEM(observations = data, covariates = covariate, FEMbasis = FEMbasis, lambda = lambda, inference.data.object=inference.data.object) # beta estimate: solution$solution$beta # tests over beta estimates p-values: solution$inference$beta$p_values # confidence intervals for beta estimates: solution$inference$beta$CI # non-parametric estimate: plot(solution$fit.FEM) # Choose lambda with GCV - stochastic grid evaluation: lambda = 10^(-2:0) solution = smooth.FEM(observations = data, covariates = covariate, FEMbasis = FEMbasis, lambda = lambda, DOF.evaluation = 'stochastic', lambda.selection.lossfunction = 'GCV') bestLambda = solution$optimization$lambda_solution # Choose lambda with GCV - Newton finite differences stochastic evaluation -: solution = smooth.FEM(observations = data, covariates = covariate, FEMbasis = FEMbasis, DOF.evaluation = 'stochastic', lambda.selection.lossfunction = 'GCV') bestLambda = solution$optimization$lambda_solution #### Smoothing with prior information about anysotropy/non-stationarity and boundary conditions #### # See Azzimonti et al. for reference to the current example data(quasicircle2D) boundary_nodes = quasicircle2D$boundary_nodes boundary_segments = quasicircle2D$boundary_segments locations = quasicircle2D$locations data = quasicircle2D$data mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) FEMbasis = create.FEM.basis(mesh) lambda = 10^-2 # Set the PDE parameters R = 2.8 K1 = 0.1 K2 = 0.2 beta = 0.5 K_func<-function(points) { output = array(0, c(2, 2, nrow(points))) for (i in 1:nrow(points)) output[,,i]=10*rbind(c(points[i,2]^2+K1*points[i,1]^2+K2*(R^2-points[i,1]^2-points[i,2]^2), (K1-1)*points[i,1]*points[i,2]), c((K1-1)*points[i,1]*points[i,2], points[i,1]^2+K1*points[i,2]^2+K2*(R^2-points[i,1]^2-points[i,2]^2))) output } b_func<-function(points) { output = array(0, c(2, nrow(points))) for (i in 1:nrow(points)) output[,i] = 10*beta*c(points[i,1],points[i,2]) output } c_func<-function(points) { rep(c(0), nrow(points)) } u_func<-function(points) { rep(c(0), nrow(points)) } PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func) # Set the boundary conditions BC = NULL BC$BC_indices = which(mesh$nodesmarkers == 1) # b.c. on the complete boundary BC$BC_values = rep(0,length(BC$BC_indices)) # homogeneus b.c. # Since the data locations are a subset of the mesh nodes for a faster solution use: dataNA = rep(NA, FEMbasis$nbasis) dataNA[mesh$nodesmarkers == 0] = data #grid evaluation solution = smooth.FEM(observations = dataNA, FEMbasis = FEMbasis, lambda = lambda, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM) # Newton's method solution = smooth.FEM(observations = dataNA, FEMbasis = FEMbasis, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM) #### Smoothing with areal data #### # See Azzimonti et al. for reference to the current exemple data(quasicircle2Dareal) incidence_matrix = quasicircle2Dareal$incidence_matrix data = quasicircle2Dareal$data mesh = quasicircle2Dareal$mesh FEMbasis = create.FEM.basis(mesh) lambda = 10^-4 # Set the PDE parameters R = 2.8 K1 = 0.1 K2 = 0.2 beta = 0.5 K_func<-function(points) { output = array(0, c(2, 2, nrow(points))) for (i in 1:nrow(points)) output[,,i]=10*rbind(c(points[i,2]^2+K1*points[i,1]^2+K2*(R^2-points[i,1]^2-points[i,2]^2), (K1-1)*points[i,1]*points[i,2]), c((K1-1)*points[i,1]*points[i,2], points[i,1]^2+K1*points[i,2]^2+K2*(R^2-points[i,1]^2-points[i,2]^2))) output } b_func<-function(points) { output = array(0, c(2, nrow(points))) for (i in 1:nrow(points)) output[,i] = 10*beta*c(points[i,1],points[i,2]) output } c_func<-function(points) { rep(c(0), nrow(points)) } u_func<-function(points) { rep(c(0), nrow(points)) } PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func) # Set the boundary conditions BC = NULL BC$BC_indices = which(mesh$nodesmarkers == 1) # b.c. on the complete boundary BC$BC_values = rep(0,length(BC$BC_indices)) # homogeneus b.c. #grid evaluation solution = smooth.FEM(observations = data, incidence_matrix = incidence_matrix, FEMbasis = FEMbasis, lambda = lambda, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM) #Newton's method solution = smooth.FEM(observations = data, incidence_matrix = incidence_matrix, FEMbasis = FEMbasis, PDE_parameters = PDE_parameters, BC = BC) plot(solution$fit.FEM) image(solution$fit.FEM)
Space-time regression with differential regularization. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.
smooth.FEM.time(locations = NULL, time_locations = NULL, observations, FEMbasis, time_mesh=NULL, covariates = NULL, PDE_parameters = NULL, BC = NULL, incidence_matrix = NULL, areal.data.avg = TRUE, FLAG_MASS = FALSE, FLAG_PARABOLIC = FALSE, FLAG_ITERATIVE = FALSE, threshold = 10^(-4), max.steps = 50, IC = NULL, search = "tree", bary.locations = NULL, family = "gaussian", mu0 = NULL, scale.param = NULL, threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15, lambda.selection.criterion = "grid", DOF.evaluation = NULL, lambda.selection.lossfunction = NULL, lambdaS = NULL, lambdaT = NULL, DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0, DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05, inference.data.object.time=NULL)smooth.FEM.time(locations = NULL, time_locations = NULL, observations, FEMbasis, time_mesh=NULL, covariates = NULL, PDE_parameters = NULL, BC = NULL, incidence_matrix = NULL, areal.data.avg = TRUE, FLAG_MASS = FALSE, FLAG_PARABOLIC = FALSE, FLAG_ITERATIVE = FALSE, threshold = 10^(-4), max.steps = 50, IC = NULL, search = "tree", bary.locations = NULL, family = "gaussian", mu0 = NULL, scale.param = NULL, threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15, lambda.selection.criterion = "grid", DOF.evaluation = NULL, lambda.selection.lossfunction = NULL, lambdaS = NULL, lambdaT = NULL, DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0, DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05, inference.data.object.time=NULL)
locations |
A matrix where each row specifies the spatial coordinates |
time_locations |
A vector containing the times of the corresponding observations in the vector |
observations |
A matrix of #locations x #time_locations with the observed data values over the spatio-temporal domain.
The spatial locations of the observations can be specified with the |
FEMbasis |
A |
time_mesh |
A vector specifying the time mesh. |
covariates |
A #observations-by-#covariates matrix where each row represents the covariates associated with the corresponding observed data value in |
PDE_parameters |
A list specifying the parameters of the PDE in the regularizing term. Default is NULL, i.e. regularization is by means of the Laplacian (stationary, isotropic case).
If the PDE is elliptic it must contain: |
BC |
A list with two vectors:
|
incidence_matrix |
A #regions-by-#triangles/tetrahedrons matrix where the element (i,j) equals 1 if the j-th triangle/tetrahedron is in the i-th region and 0 otherwise.
This is only for areal data. In case of pointwise data, this parameter is set to |
areal.data.avg |
Boolean. It involves the computation of Areal Data. If |
FLAG_MASS |
Boolean. This parameter is considered only for separable problems i.e. when |
FLAG_PARABOLIC |
Boolean. If |
FLAG_ITERATIVE |
Boolean. If |
threshold |
This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.
Default value |
max.steps |
This parameter is used to limit the maximum number of iteration.
Default value |
IC |
Initial condition needed in case of parabolic problem i.e. when |
search |
a flag to decide the search algorithm type (tree or naive or walking search algorithm). |
bary.locations |
A list with three vectors:
|
family |
This parameter specify the distribution within exponential family used for GLM model.
The following distribution are implemented: "binomial", "exponential", "gamma", "poisson", "gaussian", "invgaussian".
The default link function for binomial is |
mu0 |
This parameter is a vector that set the starting point for FPIRLS algorithm. It represent an initial guess of the location parameter.
Default is set to observation for non binary distribution while equal to |
scale.param |
Dispersion parameter of the chosen distribution. This is only required for "gamma", "gaussian", "invgaussian". User may specify the parameter as a positive real number. If the parameter is not supplied, it is estimated from data according to Wilhelm Sangalli 2016. |
threshold.FPIRLS |
This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.FPIRLS.
Default value |
max.steps.FPIRLS |
This parameter is used to limit the maximum number of iteration.
Default value |
lambda.selection.criterion |
This parameter is used to select the optimization method related to smoothing parameter |
DOF.evaluation |
This parameter is used to identify if and how degrees of freedom computation has to be performed.
The following possibilities are allowed: NULL, 'exact' and 'stochastic'
In the former case no degree of freedom is computed, while the other two methods enable computation.
Stochastic computation of DOFs may be slightly less accurate than its deterministic counterpart, but is highly suggested for meshes of more than 5000 nodes, being fairly less time consuming.
Default value |
lambda.selection.lossfunction |
This parameter is used to understand if some loss function has to be evaluated.
The following possibilities are allowed: NULL and 'GCV' (generalized cross validation)
The former case is that of |
lambdaS |
A scalar or vector of spatial smoothing parameters. |
lambdaT |
A scalar or vector of temporal smoothing parameters. |
DOF.stochastic.realizations |
This parameter is considered only when |
DOF.stochastic.seed |
This parameter is considered only when |
DOF.matrix |
Matrix of degrees of freedom. This parameter can be used if the DOF.matrix corresponding to |
GCV.inflation.factor |
Tuning parameter used for the estimation of GCV. Default value |
lambda.optimization.tolerance |
Tolerance parameter, a double between 0 and 1 that fixes how much precision is required by the optimization method: the smaller the parameter, the higher the accuracy.
Used only if |
inference.data.object.time |
An |
A list with the following variables:
fit.FEM.timeA FEM.time object that represents the fitted spatio-temporal field.
PDEmisfit.FEM.timeA FEM.time object that represents the misfit of the penalized PDE.
betaIf covariates is not NULL, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. The jth column represents the vector of regression coefficients when
the smoothing parameter is equal to lambda[j].
edfIf GCV is TRUE, a scalar or matrix with the trace of the smoothing matrix for each combination of the smoothing parameters specified in lambdaS and lambdaT.
stderrIf GCV is TRUE, a scalar or matrix with the estimate of the standard deviation of the error for each combination of the smoothing parameters specified in lambdaS and lambdaT.
GCVIf GCV is TRUE, a scalar or matrix with the value of the GCV criterion for each combination of the smoothing parameters specified in lambdaS and lambdaT.
bestlambdaIf GCV is TRUE, a 2-elements vector with the indices of smoothing parameters returning the lowest GCV
ICestimatedIf FLAG_PARABOLIC is TRUE and IC is NULL, a list containing a FEM object with the initial conditions, the value of the smoothing parameter lambda returning the lowest GCV and, in presence of covariates, the estimated beta coefficients
bary.locationsA barycenter information of the given locations if the locations are not mesh nodes.
inferenceA list set only if a well defined inferenceDataObjectTime is passed as parameter to the function; contains all inference outputs required:
p_valueslist of lists set only if at least one p-value is required; contains the p-values divided by implementation:
waldlist containing all the Wald p-values required, in the same order of the type list in inference.data.object.time.
If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of coeff matrix in inference.data.object.time.
speckmanlist containing all the Speckman p-values required, in the same order of the type list in inference.data.object.time.
If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of coeff matrix in inference.data.object.time.
eigen_sign_fliplist containing all the Eigen-Sign-Flip p-values required, in the same order of the type list in inference.data.object.time.
If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of coeff matrix in inference.data.object.time.
CIlist of lists set only if at least one confidence interval is required; contains the confidence intervals divided by implementation:
waldlist containing all the Wald confidence intervals required, in the same order of the type list in inference.data.object.time.
Each item is a matrix with 3 columns and p rows, p being the number of rows of coeff matrix in inference.data.object.time; each row is the CI for the corresponding row of coeff matrix.
speckmanlist containing all the Speckman confidence intervals required, in the same order of the type list in inference.data.object.time.
Each item is a matrix with 3 columns and p rows, p being the number of rows of coeff matrix in inference.data.object.time; each row is the CI for the corresponding row of coeff matrix.
#' @references Arnone, E., Azzimonti, L., Nobile, F., & Sangalli, L. M. (2019). Modeling spatially dependent functional data via regression with differential regularization. Journal of Multivariate Analysis, 170, 275-295. Bernardi, M. S., Sangalli, L. M., Mazza, G., & Ramsay, J. O. (2017). A penalized regression model for spatial functional data with application to the analysis of the production of waste in Venice province. Stochastic Environmental Research and Risk Assessment, 31(1), 23-38. Federico Ferraccioli, Laura M. Sangalli & Livio Finos (2022). Some first inferential tools for spatial regression with differential regularization. Journal of Multivariate Analysis, 189, 104866.
library(fdaPDE) data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations time_locations = seq(0,1,length.out = 5) mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) space_time_locations = cbind(rep(time_locations,each=nrow(mesh$nodes)), rep(mesh$nodes[,1],5),rep(mesh$nodes[,2],5)) FEMbasis = create.FEM.basis(mesh) lambdaS = 10^-1 lambdaT = 10^-1 data = fs.test(space_time_locations[,2], space_time_locations[,3])*cos(pi*space_time_locations[,1]) + rnorm(nrow(space_time_locations), sd = 0.5) data = matrix(data, nrow = nrow(mesh$nodes), ncol = length(time_locations), byrow = TRUE) solution = smooth.FEM.time(observations = data, time_locations = time_locations, FEMbasis = FEMbasis, lambdaS = lambdaS, lambdaT = lambdaT) plot(solution$fit.FEM)library(fdaPDE) data(horseshoe2D) boundary_nodes = horseshoe2D$boundary_nodes boundary_segments = horseshoe2D$boundary_segments locations = horseshoe2D$locations time_locations = seq(0,1,length.out = 5) mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments) space_time_locations = cbind(rep(time_locations,each=nrow(mesh$nodes)), rep(mesh$nodes[,1],5),rep(mesh$nodes[,2],5)) FEMbasis = create.FEM.basis(mesh) lambdaS = 10^-1 lambdaT = 10^-1 data = fs.test(space_time_locations[,2], space_time_locations[,3])*cos(pi*space_time_locations[,1]) + rnorm(nrow(space_time_locations), sd = 0.5) data = matrix(data, nrow = nrow(mesh$nodes), ncol = length(time_locations), byrow = TRUE) solution = smooth.FEM.time(observations = data, time_locations = time_locations, FEMbasis = FEMbasis, lambdaS = lambdaS, lambdaT = lambdaT) plot(solution$fit.FEM)
A dataset with information about the connectivity matrix and the nodes locations of a sphere geometry. It containes:
nodes. A #nodes-by-3 matrix specifying the locations of each node.
tetrahedrons. A #tetrahedrons-by-4 matrix specifying the indices of the nodes in each tetrahedron.
This dataset can be used to create a MESH.3D object with the function create.MESH.3D