smooth.FEM.PDE.basis: Spatial regression with differential regularization: anysotropic case (elliptic PDE)

Description

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.

Usage

smooth.FEM.PDE.basis(locations = NULL, observations, FEMbasis, 
       lambda, PDE_parameters, covariates = NULL, BC = NULL, GCV = FALSE, 
       CPP_CODE = TRUE)

Arguments

locations

A #observations-by-2 matrix where each row specifies the spatial coordinates x and y of the corresponding observations in the vector observations. This parameter can be NULL. In this case the spatial coo

observations

A vector of length #observations with the observed data values over the domain. The locations of the observations can be specified with the locations argument. Otherwise if only the vector of observations is given, these are consider to be

FEMbasis

A FEMbasis object describing the Finite Element basis, as created by create.FEM.basis.

lambda

A scalar or vector of smoothing parameters.

PDE_parameters

A list specifying the parameters of the elliptic PDE in the regularizing term: K, a 2-by-2 matrix of diffusion coefficients. This induces an anisotropic smoothing with a preferential direction that corresponds to the first eigenvector of the

covariates

A #observations-by-#covariates matrix where each row represents the covariates associated with the corresponding observed data value in observations.

A list with two vectors: BC_indices, a vector with the indices in nodes of boundary nodes where a Dirichlet Boundary Condition should be applied; BC_values, a vector with the values that the spatial field must take

GCV

Boolean. If TRUE the following quantities are computed: the trace of the smoothing matrix, the estimated error standard deviation, and the Generalized Cross Validation criterion, for each value of the smoothing parameter specified in

CPP_CODE

Boolean. If TRUE the computation relies on the C++ implementation of the algorithm. This usually ensures a much faster computation.

Value

A list with the following variables:
fit.FEMA FEM object that represents the fitted spatial field.
PDEmisfit.FEMA FEM object that represents the PDE misfit for the estimated spatial field.
betaIf covariates is not NULL, a matrix with number of rows equal to the number of covariates and numer 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 vector with the trace of the smoothing matrix for each value of the smoothing parameter specified in lambda.
stderrIf GCV is TRUE, a scalar or vector with the estimate of the standard deviation of the error for each value of the smoothing parameter specified in lambda.
GCVIf GCV is TRUE, a scalar or vector with the value of the GCV criterion for each value of the smoothing parameter specified in lambda.

References

Azzimonti, L., Sangalli, L.M., Secchi, P., Domanin, M., and Nobile, F., 2014. Blood flow velocity field estimation via spatial regression with PDE penalization Blood flow velocity field estimation via spatial regression with PDE penalization. DOI. 10.1080/01621459.2014.946036. Azzimonti, L., Nobile, F., Sangalli, L.M., and Secchi, P., 2014. Mixed Finite Elements for Spatial Regression with PDE Penalization. SIAM/ASA Journal on Uncertainty Quantification, 2(1), pp.305-335.

Examples

Run this code

# Load the mesh and plot it
data(mesh.2D.simple)
plot(mesh.2D.simple)
# Create a vector of noisy samples of an underlying spatial field, 
# located over the nodes of the mesh
observations = sin(pi*mesh.2D.simple$nodes[,1]) + rnorm(n = nrow(mesh.2D.simple$nodes), sd = 0.1)
# Create the FEM basis object
FEMbasis = create.FEM.basis(mesh.2D.simple)

# Set a vector of smoothing coefficients
lambda = c(10^-4, 1, 10^4)

# Set the anysotropic smoothing matrix K
PDE_parameters_anys = list(K = matrix(c(0.01,0,0,1), nrow = 2), b = c(0,0), c = 0)
# Estimate one field for each smoothing parameter and plot these
FEM_CPP_PDE = smooth.FEM.PDE.basis(observations = observations, 
                                   FEMbasis = FEMbasis, lambda = lambda, 
                                   PDE_parameters = PDE_parameters_anys)
plot(FEM_CPP_PDE$fit.FEM)

# Evaluate solution in three points
eval.FEM(FEM_CPP_PDE$fit.FEM, locations = rbind(c(0,0),c(0.5,0),c(-2,-2)))