Welcome to the fdaPDE 2.0’s alpha-testing!¶

library(fdaPDE2)
data("unit_square", package = "fdaPDE2")
unit_square <- Mesh(unit_square)

## the functional space of finite element functions of order 1, over the unit square
Vh <- FunctionSpace(unit_square, fe_order = 1)
basis <- Vh$get_basis() ## recover the (lagrangian) basis system of Vh

## evaluate the basis system on a given set of locations (compute Psi matrix)
locations <- unit_square$nodes
Psi <- basis$eval(type = "pointwise", locations)

## integrate a function (expressed as basis expansion on Vh) over the domain
f <- function(p) { p[,1]^2 + p[,2]^2 } ## x^2 + y^2
Vh$integrate(f)

library(fdaPDE2)
data("unit_square", package = "fdaPDE2")
unit_square <- Mesh(unit_square)

## the functional space of finite element functions of order 1, over the unit square
Vh <- FunctionSpace(unit_square, fe_order = 1)
f  <- Function(Vh) ## a generic element of the space Vh

## compose the differential operator (in strong form)
Lf <- -laplace(f) + dot(c(1,1), grad(f)) ## a costant coefficients advection-diffusion problem
## define the forcing term
u <- function(points) { return(rep(1, times = nrow(points))) }

## create your penalty
penalty <- PDE(Lf, u)

## general linear second order parabolic operator
Lf <- dt(f) - div(K * grad(f)) + dot(b, grad(f)) + c * f

library(fdaPDE2)
data("unit_square", package = "fdaPDE2")
unit_square <- Mesh(unit_square)

data_frame <- ## obtain your data in some way...
## currently, the only requirement is data_frame to be a data.frame object, e.g.,

##            y         x1         x2
## 1 -0.04044303  0.1402058 0.00000000
## 2  0.15079619  1.1989599 0.03447593
## 3  0.02391597 -2.3299685 0.06891086
## 4  0.38927632  0.5709451 0.10326387
## 5  0.39417457  2.7482761 0.13749409
## 6  0.33297548  1.7080400 0.17156085

## a model is first described in an "abstract" way, think for instance to the continuous
## functional J(f, \beta) we minimize when solving a smoothing problem

## a nonparametric spatial regression model
model <- SRPDE(y ~ f, domain = unit_square, data = data_frame, lambda = 1e-6)

## this will inject in the current environment a Function object named f, representing the
## unknown spatial field and defined on a FunctionSpace(unit_square, fe_order = 1)
## the name of the spatial field supplied in formula can be any

model$fit() ## fits the model

## we can describe a semi-parametric spatial regression model as
model <- SRPDE(y ~ x1 + f, domain = unit_square, data = data_frame, lambda = 1e-6)
model$fit()

## define the differential penalty
f <- Function(FunctionSpace(unit_square, fe_order = 1))
Lf <- -laplace(f) ## simple laplacian penalty
u <- function(points) { return(rep(0, times = nrow(points))) }

## supply the penalty to the model and fit
model <- SRPDE(y ~ f, penalty = PDE(Lf, u), data = data_frame, lambda = 1e-6)
model$fit()

model <- SRPDE(y ~ f, penalty = PDE(Lf, u), data = data_frame, lambda = 1e-6)

library(fdaPDE2)
data("unit_square", package = "fdaPDE2")
unit_square <- Mesh(unit_square)

data_frame <- ## obtain your data in some way...
model <- SRPDE(y ~ f, domain = unit_square, data = data_frame) ## do not fix any lambda...

## ... because you might want to select it via GCV minimization!
lambda_grid = 10^seq(-4, -3, by = 0.1)
model$fit(
   lambda = gcv(optimizer = "grid", lambda = lambda_grid)
)

library(fdaPDE2)
data("unit_square", package = "fdaPDE2")
unit_square <- Mesh(unit_square)

data_frame <- ## obtain your data in some way...
model <- SRPDE(y ~ f, domain = unit_square, data = data_frame, family = "poisson")

lambda_grid = 10^seq(-4, -3, by = 0.1)
model$fit(
   lambda = gcv(optimizer = "grid", lambda = lambda_grid),
   fprils_params = list( ## customize FPIRLS execution
      max_iterations = 100,
      tolerance = 1e-6
   )
)

## you are not forced to set all the parameters, if you do not set some (or all) of them,
## there will always be a default

## stochastic approximation of degrees of freedom, fast but approximate
model$fit(
   lambda = gcv(
      edf_computation = "stochastic",
      seed = 143547654,    ## set seed in internal stochastic engine (for reproducibility)
      n_mc_samples = 500,  ## number of mc samples in stochastic algorithm
      optimizer = "...",
      ... ## additional arguments forwarded to the optimizer
)

## exact approximation of degrees of freedom, slow but exact
model$fit(
   lambda = gcv(
      edf_computation = "exact",
      optimizer = "...",
      ... ## additional arguments forwarded to the optimizer
)