Nonparametric density estimation with differential regularization

Description

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.

Usage

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)

Arguments

Value

A list with the following variables:

References

• 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.

Examples

library(fdaPDE)

## Create a 2D mesh over a squared domain
Xbound <- seq(-3, 3, length.out = 10)
Ybound <- seq(-3, 3, length.out = 10)
grid_XY <- expand.grid(Xbound, Ybound)
Bounds <- grid_XY[(grid_XY$Var1 %in% c(-3, 3)) | (grid_XY$Var2 %in% c(-3, 3)), ]
mesh <- create.mesh.2D(nodes = Bounds, order = 1)
mesh <- refine.mesh.2D(mesh, maximum_area = 0.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")

Nonparametric spatio-temporal density estimation with differential regularization

Description

This function implements a nonparametric spatio-temporal density estimation method with differential regularization (given by the sum of the square of the L2 norm of the laplacian of the density function and the square of the L2 norm of the second- order time-derivative), when points are located over a planar mesh. The computation relies only on the C++ implementation of the algorithm.

Usage

DE.FEM.time(data, data_time, FEMbasis, mesh_time, lambda, lambda_time, scaling=NULL,
                   fvec=NULL, heatStep=0.1, heatIter=10, stepProposals=NULL, tol1=1e-4,
                   tol2=0, print=FALSE, nfolds=0, nsimulations=500,
                   step_method="Fixed_Step", direction_method="BFGS",
                   preprocess_method="NoCrossValidation", search="tree",
                   isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE,
                   inference = FALSE)

Arguments

Value

A list with the following variables:

Examples

library(fdaPDE)

## Create a 2D mesh over a squared domain
Xbound <- seq(-3, 3, length.out = 10)
Ybound <- seq(-3, 3, length.out = 10)
grid_XY <- expand.grid(Xbound, Ybound)
Bounds <- grid_XY[(grid_XY$Var1 %in% c(-3, 3)) | (grid_XY$Var2 %in% c(-3, 3)), ]
mesh <- create.mesh.2D(nodes = Bounds, order = 1)
mesh <- refine.mesh.2D(mesh, maximum_area = 0.25, minimum_angle = 20)
FEMbasis <- create.FEM.basis(mesh)

## Create a 1D time mesh over a (non-negative) interval
mesh_time <- seq(0, 1, length.out=3)

## Generate data
n <- 50
set.seed(10)
x <- rnorm(n,0,2)
y <- rnorm(n,0,2)
locations <- cbind(x,y)
times <- runif(n,0,1)
data <- cbind(locations, times)

t <- 0.5 # time instant in which to evaluate the solution

plot(mesh)
sample <- data[abs(data[,3]-t)<0.05,1:2]
points(sample, col="red", pch=19, cex=1, main=paste('Sample | ', t-0.05,' < t < ', t+0.05))

## Density Estimation
lambda <- 0.1
lambda_time <- 0.001
sol <- DE.FEM.time(data = locations, data_time = times, FEMbasis = FEMbasis, mesh_time = mesh_time,
                   lambda = lambda, lambda_time = lambda_time, nsimulations=300, inference=TRUE)

## Visualization
n = 100
X <- seq(-3, 3, length.out = n)
Y <- seq(-3, 3, length.out = n)
grid <- expand.grid(X, Y)

FEMfunction = FEM.time(sol$g, mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation <- eval.FEM.time(FEM.time = FEMfunction, locations = grid, time.instants = t)
FEMfunction_L = FEM.time(sol$g_CI_L, mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation_L <- eval.FEM.time(FEM.time = FEMfunction_L, locations = grid, time.instants = t)
FEMfunction_U = FEM.time(sol$g_CI_U, mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation_U <- eval.FEM.time(FEM.time = FEMfunction_U, locations = grid, time.instants = t)

image2D(x = X, y = Y, z = matrix(exp(evaluation_L), n, n), col = heat.colors(100),
        xlab = "x", ylab = "y", contour = list(drawlabels = FALSE),
        main = paste("Estimated CI lower bound at t = ", t), zlim=c(0,0.3), asp = 1)
image2D(x = X, y = Y, z = matrix(exp(evaluation), n, n), col = heat.colors(100),
        xlab = "x", ylab = "y", contour = list(drawlabels = FALSE),
        main = paste("Estimated density at t = ", t), zlim=c(0,0.3), asp = 1)
image2D(x = X, y = Y, z = matrix(exp(evaluation_U), n, n), col = heat.colors(100),
        xlab = "x", ylab = "y", contour = list(drawlabels = FALSE),
        main = paste("Estimated CI upper bound at t = ", t), zlim=c(0,0.3), asp = 1)

Density initialization

Description

This function implements two methods for the density initialization procedure.

Usage

DE.heat.FEM(data, FEMbasis, lambda=NULL, heatStep=0.1, heatIter=500, 
            init="Heat", nFolds=5, search = "tree")

Arguments

Value

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.

Examples

library(fdaPDE)

## Create a 2D mesh over a squared domain
Xbound <- seq(-3, 3, length.out = 10)
Ybound <- seq(-3, 3, length.out = 10)
grid_XY <- expand.grid(Xbound, Ybound)
Bounds <- grid_XY[(grid_XY$Var1 %in% c(-3, 3)) | (grid_XY$Var2 %in% c(-3, 3)), ]
mesh <- create.mesh.2D(nodes = Bounds, order = 1)
mesh <- refine.mesh.2D(mesh, maximum_area = 0.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))

Spatio-temporal density initialization

Description

This function implements two methods for the density initialization procedure.

Usage

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)

Arguments

Value

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.

Examples

library(fdaPDE)

## Create a 2D mesh over a squared domain
Xbound <- seq(-3, 3, length.out = 10)
Ybound <- seq(-3, 3, length.out = 10)
grid_XY <- expand.grid(Xbound, Ybound)
Bounds <- grid_XY[(grid_XY$Var1 %in% c(-3, 3)) | (grid_XY$Var2 %in% c(-3, 3)), ]
mesh <- create.mesh.2D(nodes = Bounds, order = 1)
mesh <- refine.mesh.2D(mesh, maximum_area = 0.25, minimum_angle = 20)
FEMbasis <- create.FEM.basis(mesh)

## Create a 1D time mesh over a (non-negative) interval
mesh_time <- seq(0, 1, length.out=3)

## Generate data
n <- 50
set.seed(10)
x <- rnorm(n,0,2)
y <- rnorm(n,0,2)
locations <- cbind(x,y)
times <- runif(n,0,1)
data <- cbind(locations, times)

t <- 0.5 # time instant in which to evaluate the solution

plot(mesh)
sample <- data[abs(data[,3]-t)<0.05,1:2]
points(sample, col="red", pch=19, cex=1, main=paste('Sample | ', t-0.05,' < t < ', t+0.05))

## Density 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)

Define a surface or spatial field by a Finite Element basis expansion

Description

This function defines a FEM object.

Usage

FEM(coeff,FEMbasis)

Arguments

Value

An FEM object. This contains a list with components coeff and FEMbasis.

Examples

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)

Define a spatio-temporal field by a Finite Element basis expansion

Description

This function defines a FEM.time object.

Usage

FEM.time(coeff,time_mesh,FEMbasis,FLAG_PARABOLIC=FALSE)

Arguments

Value

A FEM.time object. This contains a list with components coeff, mesh_time, FEMbasis and FLAG_PARABOLIC.

Examples

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)

Smooth Functional Principal Component Analysis

Description

This function implements a smooth functional principal component analysis over a planar mesh, a smooth manifold or a volume.

Usage

FPCA.FEM(locations = NULL, datamatrix, FEMbasis, lambda, nPC = 1, validation = NULL,
                NFolds = 5,GCVmethod = "Stochastic", nrealizations = 100, search = "tree",
                bary.locations = NULL)

Arguments

Value

A list with the following variables:

loadings.FEM

A FEM object that represents the L^2-normalized functional loadings for each Principal Component computed.

scores

A #samples-by-#PrincipalComponents matrix that represents the unnormalized scores or PC vectors.

lambda

A vector of length #PrincipalComponents with the values of the smoothing parameter lambda chosen for that Principal Component.

variance_explained

A vector of length #PrincipalComponents where each value represent the variance explained by that component.

cumsum_percentage

A vector of length #PrincipalComponents containing the cumulative percentage of the variance explained by the first components.

bary.locations

A barycenter information of the given locations if the locations are not mesh nodes.

References

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.

Examples

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)

Covariate test function for the horseshoe domain

Description

Implements a finite area test function the horseshoe domain.

Usage

covs.test(x, y)

Arguments

Value

Returns function evaluations.

Create a FEM basis

Description

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.

Usage

create.FEM.basis(mesh, saveTree = FALSE)

Arguments

Value

A FEMbasis object. This contains the mesh, along with some additional quantities:

order

Either "1" or "2" for the 2D and 2.5D case, and "1" for the 3D case. Order of the Finite Element basis.

nbasis

Scalar. The number of basis.

See Also

create.mesh.2D, create.mesh.2.5D,create.mesh.3D

Examples

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

Description

Create a 1.5D linear network mesh

Usage

create.mesh.1.5D(nodes, edges = NULL, order = 1, nodesattributes = NULL)

Arguments

Value

An object of the class mesh.1.5D 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

A matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.

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 nodes, indicating the nodes where the edge starts from and ends to.

neighbors

A #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.

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.

Create a mesh.2.5D object from the nodes locations and the connectivity matrix

Description

Create a mesh.2.5D object from the nodes locations and the connectivity matrix

Usage

create.mesh.2.5D(
  nodes,
  triangles = NULL,
  order = 1,
  nodesattributes = NULL,
  segments = NULL,
  holes = NULL
)

Arguments

Value

An object of the class mesh.2.5D with the following output:

nodes

A #nodes-by-3 matrix containing the x, y, z 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

A matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.

triangles

A #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.

segmentsmarkers

A 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.

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 nodes, indicating 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 edge is a boundary edge; an entry '0' indicates that the corresponding edge is not a boundary edge.

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-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.

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.

Examples

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)

Create a 2D triangular mesh

Description

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.

Usage

create.mesh.2D(nodes, nodesattributes = NA, segments = NA, holes = NA,
                     triangles = NA, order = 1, verbosity = 0)

Arguments

Value

An object of the class mesh.2D 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

A matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.

triangles

A #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.

segmentsmarkers

A 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.

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 nodes, indicating 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 edge is a boundary edge; an entry '0' indicates that the corresponding edge is not a boundary edge.

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.

See Also

refine.mesh.2D, create.FEM.basis

Examples

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.

Create a mesh.3D object from the connectivity matrix and nodes locations

Description

Create a mesh.3D object from the connectivity matrix and nodes locations

Usage

create.mesh.3D(
  nodes,
  tetrahedrons,
  order = 1,
  nodesattributes = NULL,
  segments = NULL,
  holes = NULL
)

Arguments

Value

An object of the class mesh.3D with the following output:

nodes

A #nodes-by-3 matrix containing the x, y, z 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

A matrix with #nodes rows containing nodes' attributes. These are passed unchanged from the input.

tetrahedrons

A #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.

segmentsmarkers

A 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.

faces

A #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.

facesmarkers

A 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.

neighbors

A #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.

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. These are passed unchanged from the input.

order

Either '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.

Examples

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)

Evaluate a FEM object at a set of point locations

Description

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.

Usage

eval.FEM(FEM, locations = NULL, incidence_matrix = NULL, search = "tree", 
                bary.locations = NULL)

Arguments

Value

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.

References

• 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.

Examples

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)

Evaluate a FEM.time object at a set of point locations

Description

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).

Usage

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)

Arguments

Value

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.

References

Devillers, O. et al. 2001. Walking in a Triangulation, Proceedings of the Seventeenth Annual Symposium on Computational Geometry

Examples

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)

Deprecated Functions

Description

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.

Usage

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, ...)

Arguments

Details

These functions are Deprecated in this release of fdaPDE, they will be marked as Defunct and removed in a future version.

Value

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:

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.FEM

A FEM object that represents the fitted spatial field.

PDEmisfit.FEM

A FEM object that represents the Laplacian of the estimated spatial field.

solution

A 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

f

Matrix, estimate of function f, first half of solution vector

g

Matrix, second half of solution vector

z_hat

Matrix, prediction of the output in the locations

beta

If 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

rmse

Estimate of the root mean square error in the locations

estimated_sd

Estiimate of the standard deviation of the error

optimization

A detailed list of optimization related data:

lambda_solution

numerical value of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"

lambda_position

integer, postion in lambda_vector of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"

GCV

numeric value of GCV in correspondence of the optimum

optimization_details

list containing further information about the optimization method used and the nature of its termination, eventual number of iterations

dof

numeric vector, value of dof for all the penalizations it has been computed, empty if not computed

lambda_vector

numeric value of the penalization factors passed by the user or found in the iterations of the optimization method

GCV_vector

numeric vector, value of GCV for all the penalizations it has been computed

time

Duration of the entire optimization computation

bary.locations

A barycenter information of the given locations if the locations are not mesh nodes.

A list with the following variables:

fit.FEM

A FEM object that represents the fitted spatial field.

PDEmisfit.FEM

A FEM object that represents the Laplacian of the estimated spatial field.

solution

A 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

f

Matrix, estimate of function f, first half of solution vector

g

Matrix, second half of solution vector

z_hat

Matrix, prediction of the output in the locations

beta

If 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

rmse

Estimate of the root mean square error in the locations

estimated_sd

Estiimate of the standard deviation of the error

optimization

A detailed list of optimization related data:

lambda_solution

numerical value of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"

lambda_position

integer, postion in lambda_vector of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"

GCV

numeric value of GCV in correspondence of the optimum

optimization_details

list containing further information about the optimization method used and the nature of its termination, eventual number of iterations

dof

numeric vector, value of dof for all the penalizations it has been computed, empty if not computed

lambda_vector

numeric value of the penalization factors passed by the user or found in the iterations of the optimization method

GCV_vector

numeric vector, value of GCV for all the penalizations it has been computed

time

Duration of the entire optimization computation

bary.locations

A barycenter information of the given locations if the locations are not mesh nodes.

A list with the following variables:

fit.FEM

A FEM object that represents the fitted spatial field.

PDEmisfit.FEM

A FEM object that represents the Laplacian of the estimated spatial field.

solution

A 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

f

Matrix, estimate of function f, first half of solution vector

g

Matrix, second half of solution vector

z_hat

Matrix, prediction of the output in the locations

beta

If 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

rmse

Estimate of the root mean square error in the locations

estimated_sd

Estiimate of the standard deviation of the error

optimization

A detailed list of optimization related data:

lambda_solution

numerical value of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"

lambda_position

integer, postion in lambda_vector of best lambda acording to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction="unused"

GCV

numeric value of GCV in correspondence of the optimum

optimization_details

list containing further information about the optimization method used and the nature of its termination, eventual number of iterations

dof

numeric vector, value of dof for all the penalizations it has been computed, empty if not computed

lambda_vector

numeric value of the penalization factors passed by the user or found in the iterations of the optimization method

GCV_vector

numeric vector, value of GCV for all the penalizations it has been computed

time

Duration of the entire optimization computation

bary.locations

A barycenter information of the given locations if the locations are not mesh nodes.

An object of the class MESH2D with the following output:

A MESH2D object representing the refined triangular mesh, with the following output:

See Also

refine.MESH.2D, create.FEM.basis

create.MESH.2D, create.FEM.basis

FELSPLINE test function

Description

Implements a finite area test function based on one proposed by Tim Ramsay (2002) proposed by Simon Wood (2008).

Usage

fs.test(x, y, r0 = 0.1, r = 0.5, l = 3, b = 1, exclude = FALSE)

Arguments

Value

Returns function evaluations, or NAs for points outside the horseshoe domain.

References

• 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.

Examples

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)

FELSPLINE 3D test function

Description

Implements a finite area test function based on one proposed by Tim Ramsay (2002) and by Simon Wood (2008) in 3D.

Usage

fs.test.3D(x, y, z, r0 = 0.25, r = 1.25, l = 5, b = 1, exclude = FALSE)

Arguments

Value

Returns function evaluations, or NAs for points outside the horseshoe domain.

Examples

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))

Horseshoe2.5D domain

Description

A mesh2.5D object with nodes and connectivity matrix of a triangular mesh of the horseshoe 2.5D domain.

Horseshoe domain

Description

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.

Hub domain

Description

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 2D FEM object

Description

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.

Usage

## S3 method for class 'FEM'
image(x, num_refinements, ...)

Arguments

Value

No return value

See Also

FEM plot.FEM

Examples

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 2D FEM.time object at a given time

Description

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.

Usage

## S3 method for class 'FEM.time'
image(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)

Arguments

Value

No return value

See Also

FEM.time image.FEM.time

Examples

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)

Class for inference data

Description

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.

Details

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.

Slots

test

A 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.

interval

A 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.

type

A 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.

component

A 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.

exact

An 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).

dim

Dimension 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_cov

Number of covariates taken into account in the linear part of the regression problem.

locations

A 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_indices

A 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_nodes

An 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;

coeff

A 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.

beta0

Vector of null hypothesis values for the linear parameters of the model. Used only if test is not 0 and component is not 2.

f0

Function representing the expression of the nonparametric component f under the null hypothesis. Used only if component is not 1.

f0_eval

Vector of f0 evaluations at the chosen test locations. It will be eventually set later in checkInferenceParameters, if nonparametric inference is required.

f_var

An 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.

quantile

Vector of quantiles needed for confidence intervals, used only if interval is not 0.

alpha

1 minus confidence level vector of sign-flipping approaches confidence intervals. Used only if interval is not 0.

n_flip

An integer representing the number of sign-flips in the case of sign-flipping approaches.

tol_fspai

A real number greater than 0 specifying the tolerance for FSPAI algorithm, in case of non-exact inference (not implemented in this version).

definition

An 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.

Constructor for inferenceDataObject class

Description

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.

Usage

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)

Arguments

Value

The output is a well defined inferenceDataObject, that can be used as input parameter in the smooth.FEM function.

Examples

obj<-inferenceDataObjectBuilder(test = 'oat', dim = 2, beta0 = rep(1,4), n_cov = 4);
obj2<-inferenceDataObjectBuilder(test = 'sim', dim = 3, component = 'nonparametric', n_cov = 3);

Class for inference data in ST case

Description

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.

Details

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.

Slots

test

A 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.

interval

A 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.

type

A vector of integers taking value 1, 2, 3 or 4 corresponding to Wald, Speckman, Eigen-Sign-Flip, Enhanced-Eigen-Sign-Flip inferential approach.

component

A 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.

exact

An 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).

dim

Dimension 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_cov

Number of covariates taken into account in the linear part of the regression problem.

locations

A 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_indices

A 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_nodes

An 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_locations

A 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.

coeff

A 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.

beta0

Vector of null hypothesis values for the linear parameters of the model. Used only if test is not 0 and component is not 2.

f0

Function representing the expression of the nonparametric component f under the null hypothesis. Used only if component is not 1.

f0_eval

Matrix of f0 evaluations at the chosen space and time locations. It will be eventually set later in checkInferenceParametersTime, if nonparametric inference is required.

f_var

An 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.

quantile

Vector of quantiles needed for confidence intervals, used only if interval is not 0.

alpha

1 minus confidence level vector of sign-flipping approaches confidence intervals. Used only if interval is not 0.

n_flip

An integer representing the number of sign-flips in the case of sign-flipping approaches.

tol_fspai

A real number greater than 0 specifying the tolerance for FSPAI algorithm, in case of non-exact inference (not implemented in this version).

definition

An 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.

Constructor for inferenceDataObjectTime class

Description

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.

Usage

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)

Arguments

Value

The output is a well defined inferenceDataObjectTime, that can be used as input parameter in the smooth.FEM function.

Examples

obj<-inferenceDataObjectTimeBuilder(test = 'oat', dim = 2, beta0 = rep(1,4), n_cov = 4);
obj2<-inferenceDataObjectTimeBuilder(test = 'sim', dim = 3, component = 'nonparametric', n_cov = 3);

Plot a FEM object

Description

Three-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).

Usage

## S3 method for class 'FEM'
plot(x, colormap = "heat.colors", num_refinements = NULL, ...)

Arguments

Value

No return value

See Also

FEM, image.FEM

Examples

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)

Plot a FEM.time object at a given time

Description

Plot 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).

Usage

## 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", ...)

Arguments

Value

No return value

See Also

FEM.time, image.FEM.time

Examples

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

Description

Plot a mesh.1.5D object, generated by create.mesh.1.5D or refine.mesh.1.5D.

Usage

## S3 method for class 'mesh.1.5D'
plot(x, ...)

Arguments

Value

No return value

Plot a mesh.2.5D object

Description

Plot the triangulation of a mesh.2.5D object, generated by create.mesh.2.5D

Usage

## S3 method for class 'mesh.2.5D'
plot(x, ...)

Arguments

Value

No return value

Examples

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

Description

Plot a mesh.2D object, generated by create.mesh.2D or refine.mesh.2D.

Usage

## S3 method for class 'mesh.2D'
plot(x, ...)

Arguments

Value

No return value

Examples

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

Description

Plot a mesh.3D object, generated by create.mesh.3D.

Usage

## S3 method for class 'mesh.3D'
plot(x, ...)

Arguments

Value

No return value

Examples

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)

Project 2D points onto 1.5D linear network mesh

Description

This function projects any 2D points onto 1.5D linear network mesh.

Usage

projection.points.1.5D(mesh, locations)

Arguments

Value

2D points projected onto 1.5D linear network mesh.

Examples

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)

Project 3D points onto 2D 2.5D triangular mesh

Description

This function projects any 3D points onto 2.5D triangular mesh.

Usage

projection.points.2.5D(mesh, locations)

Arguments

Value

3D points projected onto 2.5D triangluar mesh.

Examples

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)

Quasicircle2D domain

Description

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.

Quasicircle2Dareal domain

Description

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.

References

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.

Create a mesh.1.5D object by splitting each edge of a given mesh into two subedges.

Description

Create a mesh.1.5D object by splitting each edge of a given mesh into two subedges.

Usage

refine.by.splitting.mesh.1.5D(mesh = NULL)

Arguments

Value

An object of class mesh.1.5D with splitted edges

Create a mesh.2.5D object by splitting each triangle of a given mesh into four subtriangles.

Description

Create a mesh.2.5D object by splitting each triangle of a given mesh into four subtriangles.

Usage

refine.by.splitting.mesh.2.5D(mesh = NULL)

Arguments

Value

An object of class mesh.2.5D with splitted triangles

Create a mesh.2D object by splitting each triangle of a given mesh into four subtriangles.

Description

Create a mesh.2D object by splitting each triangle of a given mesh into four subtriangles.

Usage

refine.by.splitting.mesh.2D(mesh = NULL)

Arguments

Value

An object of class mesh.2D with splitted triangles

Create a mesh.3D object by splitting each tetrahedron of a given mesh into eight subtetrahedrons.

Description

Create a mesh.3D object by splitting each tetrahedron of a given mesh into eight subtetrahedrons.

Usage

refine.by.splitting.mesh.3D(mesh = NULL)

Arguments

Value

An object of class mesh.3D with splitted tetrahedrons

Refine 1.5D mesh

Description

Refine 1.5D mesh

Usage

refine.mesh.1.5D(mesh, delta)

Arguments

Value

An object of class mesh.1.5D with refined edges

Refine a 2D triangular mesh

Description

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.

Usage

refine.mesh.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity)

Arguments

Value

A mesh.2D 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 order = 1) or #triangles-by-6 (when order = 2) matrix.

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 edge is a boundary edge; an entry '0' indicates that the corresponding edge is not a boundary edge.

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.

See Also

create.mesh.2D, create.FEM.basis

Examples

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)

Spatial regression with differential regularization

Description

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

Usage

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)

Arguments

Value

A list with the following variables:

fit.FEM

A FEM object that represents the fitted spatial field.

PDEmisfit.FEM

A FEM object that represents the Laplacian of the estimated spatial field.

solution

A 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.

f

Matrix, estimate of function f, first half of solution vector.

g

Matrix, second half of solution vector.

z_hat

Matrix, prediction of the output in the locations.

beta

If 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.

rmse

Estimate of the root mean square error in the locations.

estimated_sd

Estimate of the standard deviation of the error.

optimization

A detailed list of optimization related data:

lambda_solution

numerical value of best lambda according to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction=NULL.

lambda_position

integer, position in lambda_vector of best lambda according to lambda.selection.lossfunction, -1 if lambda.selection.lossfunction=NULL.

GCV

numeric value of GCV in correspondence of the optimum.

optimization_details

list containing further information about the optimization method used and the nature of its termination, eventual number of iterations.

dof

vector of positive numbers, DOFs for all the lambdas in lambda_vector, empty or invalid if not computed.

lambda_vector

vector of positive numbers: penalizations either passed by the user or found in the iterations of the optimization method.

GCV_vector

vector of positive numbers, GCV values for all the lambdas in lambda_vector

time

Duration of the entire optimization computation.

bary.locations

Barycenter information of the given locations, if the locations are not mesh nodes.

GAM_output

A list of GAM related data:

fn_hat

A 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_minima

A vector of the same length of lambda, containing the reached minima for each value of the smoothing parameter.

variance.est

A vector which returns the variance estimates for the Generative Additive Models.

inference

A list set only if a well defined inferenceDataObject is passed as parameter to the function; contains all inference outputs required:

p_values

list of lists set only if at least one p-value is required; contains the p-values divided by implementation:

wald

list 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.

speckman

list 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_flip

list 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.

CI

list of lists set only if at least one confidence interval is required; contains the confidence intervals divided by implementation:

wald

list 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.

speckman

list 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.

References

• 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.

Examples

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

Description

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.

Usage

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)

Arguments

Value

A list with the following variables:

fit.FEM.time

A FEM.time object that represents the fitted spatio-temporal field.

PDEmisfit.FEM.time

A FEM.time object that represents the misfit of the penalized PDE.

beta

If 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].

edf

If 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.

stderr

If 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.

GCV

If 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.

bestlambda

If GCV is TRUE, a 2-elements vector with the indices of smoothing parameters returning the lowest GCV

ICestimated

If 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.locations

A barycenter information of the given locations if the locations are not mesh nodes.

inference

A list set only if a well defined inferenceDataObjectTime is passed as parameter to the function; contains all inference outputs required:

p_values

list of lists set only if at least one p-value is required; contains the p-values divided by implementation:

wald

list 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.

speckman

list 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_flip

list 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.

CI

list of lists set only if at least one confidence interval is required; contains the confidence intervals divided by implementation:

wald

list 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.

speckman

list 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

#' @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.

Examples

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)

Sphere3Ddata

Description

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