| Title: | Spatio-Temporal (Hero) Sandwich Smoother |
|---|---|
| Description: | An implementation of the sandwich smoother proposed in Fast Bivariate Penalized Splines by Xiao et al. (2012) <doi:10.1111/rssb.12007>. A hero is a specific type of sandwich. Dictionary.com (2018) <https://www.dictionary.com> describes a hero as: a large sandwich, usually consisting of a small loaf of bread or long roll cut in half lengthwise and containing a variety of ingredients, as meat, cheese, lettuce, and tomatoes. Also implements the spatio-temporal sandwich smoother of French and Kokoszka (2021) <doi:10.1016/j.spasta.2020.100413>. |
| Authors: | Joshua French |
| Maintainer: | Joshua French <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.6 |
| Built: | 2024-11-04 05:21:31 UTC |
| Source: | https://github.com/cran/hero |
adjacent attempts to find the point(s) adjacent
(closest) to each point. The data are implicitly assumed
to be on a grid, otherwise this function isn't very
useful. Distances between each point and other points in
coords are computed and then rounded using the
round function. Let k denote
the minimum distance between a reference point and all
other points. A point is adjacent to the reference point
if (after rounding), it's distance from the reference
point matches the minimum distance k.
adjacent(coords, longlat = FALSE, digits = 1)adjacent(coords, longlat = FALSE, digits = 1)
coords |
A two-dimensional matrix-like object with non-NULL dimensions. |
longlat |
A logical value indicating whether Great
Circle distances should be used ( |
digits |
The number of digits to use when applying
|
digits is the number of digits used by
round in the rounding process.
A hero_adjacent object. This is simply a
list with elements nbrs and coords.
nbrs is a list specifying the adjacent points
for each point. coords is simply the original
coords supplied to the function and is retained
for plotting purposes.
# basic coordinates coords = expand.grid(1:4, 1:4) # plot coordinates to see relationships plot(coords, type = "n") text(coords) a = adjacent(coords, digits = 1) plot(a)# basic coordinates coords = expand.grid(1:4, 1:4) # plot coordinates to see relationships plot(coords, type = "n") text(coords) a = adjacent(coords, digits = 1) plot(a)
starray
Convert a three-dimensional spatio-temporal array into
an starray object. The first two dimensions are
assumed to relate to gridded spatial positions.
as.starray(x) starray(x) as_starray(x)as.starray(x) starray(x) as_starray(x)
x |
A three-dimensional array |
An starray object.
star = as.starray(tasmax) class(star)star = as.starray(tasmax) class(star)
sts classConvert a numeric three-dimensional array or
two-dimensional matrix-like object to an
sts (spatial time series) object. If x is
a three-dimensional array, the first two dimensions are
assumed to relate to gridded spatial positions. If
x has only two dimensions, each row is a time
series for a specific location. Each column is a
realization of a geostatistical process at a specific
time.
as.sts(x) sts(x) as_sts(x)as.sts(x) sts(x) as_sts(x)
x |
A matrix-like object with 2 dimensions or an array with 3 dimensions. |
This method has been tested with objects of class
matrix, data.frame,
array, and
Matrix-class. It should be
possible for x to have a different class as long
as the object has a loaded as.matrix
method, which is used in this function.
An sts object.
# 3d array to sts sts = as.sts(tasmax) class(sts) # extract a subset of tasmax to produce an sts x = matrix(c(tasmax[50:60, 50:60, ]), ncol = 30) sts = as.sts(x) class(sts) sts = as.sts(as.array(x)) class(sts) sts = as.sts(Matrix::Matrix(x)) class(sts) sts = as.sts(as.data.frame(x)) class(sts)# 3d array to sts sts = as.sts(tasmax) class(sts) # extract a subset of tasmax to produce an sts x = matrix(c(tasmax[50:60, 50:60, ]), ncol = 30) sts = as.sts(x) class(sts) sts = as.sts(as.array(x)) class(sts) sts = as.sts(Matrix::Matrix(x)) class(sts) sts = as.sts(as.data.frame(x)) class(sts)
Assemble computations from a spline-related object
x in order to implement the sandwich smoother.
This is essentially an internal function, but could be
useful for performing a manual implementation of the
sandwich smoother.
assemble(object, ...) ## S3 method for class 'hero_bspline' assemble(object, x, m = 2, sparse = TRUE, ...) ## S3 method for class 'hero_radspline' assemble(object, x, m = 2, sparse = TRUE, spdiffpen = TRUE, digits = 1, ...) ## S3 method for class 'list' assemble(object, x, m = 2, sparse = TRUE, spdiffpen = TRUE, digits = 1, ...)assemble(object, ...) ## S3 method for class 'hero_bspline' assemble(object, x, m = 2, sparse = TRUE, ...) ## S3 method for class 'hero_radspline' assemble(object, x, m = 2, sparse = TRUE, spdiffpen = TRUE, digits = 1, ...) ## S3 method for class 'list' assemble(object, x, m = 2, sparse = TRUE, spdiffpen = TRUE, digits = 1, ...)
object |
A spline-related object (e.g, a
|
... |
Not implemented |
x |
Values at which to evaluate the basis functions.
This
should be a numeric vector if |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
spdiffpen |
A logical value indicating whether
|
digits |
The number of digits to use when applying
|
A list with the necessary components (ingredients)
# construct b-spline object1 = bspline(nbasis = 10) # sequence to evaluate x1 = seq(0, 1, len = 11) # assemble b-spline information spline1 = assemble(object1, x1) # assemble radial spline information border = border.grid(lon, lat) object2 = radspline(nknots = 16, border) x2 = cbind(c(lon), c(lat)) spline2 = assemble(object2, x = x2) # assemble for list of splines object = list(object1, object2) x = list(x1, x2) splines = assemble(object, x)# construct b-spline object1 = bspline(nbasis = 10) # sequence to evaluate x1 = seq(0, 1, len = 11) # assemble b-spline information spline1 = assemble(object1, x1) # assemble radial spline information border = border.grid(lon, lat) object2 = radspline(nknots = 16, border) x2 = cbind(c(lon), c(lat)) spline2 = assemble(object2, x = x2) # assemble for list of splines object = list(object1, object2) x = list(x1, x2) splines = assemble(object, x)
border.grid determines the border for
data on a grid. x and y must define a
regular or irregular grid. See Details.
border.grid(x, y, proj4string) border_grid(x, y, proj4string) borderGrid(x, y, proj4string) BorderGrid(x, y, proj4string)border.grid(x, y, proj4string) border_grid(x, y, proj4string) borderGrid(x, y, proj4string) BorderGrid(x, y, proj4string)
x |
A vector or matrix of x coordinates. See Details. |
y |
A vector or matrix of y coordinates. See Details. |
proj4string |
A projection string of class
|
A regular grid is defined by ascending numeric vectors
x and y. A vector x is ascending if
x[i] < x[j] for i < j.
An irregular grid is defind by ascending matrices.
A matrix x is ascending if x[i, j] < x[i, l]
for j < l and if x[i, j] < x[k, j]
and j < k.
A SpatialPolygons object.
Joshua French
# create x and y defining square border x = seq(min(lon), max(lon), length = 60) y = seq(min(lat), max(lat), length = 80) border = border.grid(x, y) sp::plot(border) # use lon and lat to define border of an irregular grid border2 = border.grid(lon, lat) sp::plot(border2)# create x and y defining square border x = seq(min(lon), max(lon), length = 60) y = seq(min(lat), max(lat), length = 80) border = border.grid(x, y) sp::plot(border) # use lon and lat to define border of an irregular grid border2 = border.grid(lon, lat) sp::plot(border2)
bspline helps define the parameters necessary for
constructing a B-spline but doesn't evaluate it.
bspline(rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, knots)bspline(rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, knots)
rangeval |
A numeric vector of length 2 defining the
interval over which the functional data object can be
evaulated. The default value is |
nbasis |
An integer specifying the number of
basis functions to construct. This is closely linked to
the number of knots ( |
nknots |
The number of *interior* knots. See Details. |
norder |
An integer specifying the order of the B-splines, which is one higher than their degree. The default is 4, which corresponds to cubic splines. |
extend |
Should the knots stop at the endpoints specified by |
knots |
A numeric vector with all knots (interior and exterior), including potentially replicated endpoints. See Details. |
The knots are assumed to be equidistant and non-repeating (except possibly at the endpoints).
The number of knots (nknots) and the number of
basis function (nbasis) are linked by the relation
nknots = nbasis - norder.
If extend = TRUE, the interior knots are
augmented by replicating the rangeval endpoints
norder times. Otherwise, the interior knots are
augmented by norder knots at each end
corresponding to the spacing of the interior knots.
The knot placement mimics the behavior of
create.bspline.basis when extend
= FALSE. Note that the number of breaks specified by
breaks in create.bspline.basis
corresponds to the number of interior knots plus 2 (the
interior knots plus the two endpoints).
If knots is specified, the function does minimial
argument checking. This is provided (mostly) for
testing purposes, though it can be used by individuals
who want more customizability of knots locations than
the equidistant spacing provided by default.
An object of class hero_bspline. It is a
list specifying the necessary B-spline parameters.
Joshua French
bspline(nbasis = 10)bspline(nbasis = 10)
The first n values of x are circulated
from the front of x to the back of x.
circulate(x, n = 1)circulate(x, n = 1)
x |
vector of values |
n |
The number of values to circulate |
The circulated vector
Joshua French
circulate(1:10, n = 2) circulate(as.list(1:10), n = 2)circulate(1:10, n = 2) circulate(as.list(1:10), n = 2)
hero_radsplines
connect joins multiple hero_radspline
objects into a single hero_radspline.
The e
connect(...)connect(...)
... |
A sequence of |
A combined hero_radspline
border = border.grid(lon, lat) s1 = radspline(nknots = 36, border = border) plot(s1) s2 = radspline(nknots = 36 * 4, border = border, width = 6) plot(s2) par(mfrow = c(1, 2)) plot(s1) plot(s2) par(mfrow = c(1, 1)) s = connect(s1, s2) plot(s)border = border.grid(lon, lat) s1 = radspline(nknots = 36, border = border) plot(s1) s2 = radspline(nknots = 36 * 4, border = border, width = 6) plot(s2) par(mfrow = c(1, 2)) plot(s1) plot(s2) par(mfrow = c(1, 1)) s = connect(s1, s2) plot(s)
prepared_list
create.prepared_list creates a
prepared_list manually. Typically, one would
simply use the prepare.list, but there are
situations where the data argument would be too
large to read into memory.
This function assumes that
the user has used the assemble function to
construct a list of the relevant assembled_splines
and manually computed Ytilde for a number of
relevant data observations and stored them in a
list. The user should also manually compute the sum of
the squared data for each data observation.
The user must also specify the dimensions of each data
set (which are assumed to be the same) as a vector and
provide the relevant set of values at which each
data object is observed. See Examples.
create.prepared_list(assembled, x, Ytilde, sum_ysq, n)create.prepared_list(assembled, x, Ytilde, sum_ysq, n)
assembled |
A list of |
x |
The list of arguments at which to evaluate each
of the splines used to construct |
Ytilde |
A list of |
sum_ysq |
A vector with the sum of squared |
n |
The dimensions of the |
A prepared list.
# generate and prepare 3d data set.seed(9) dat = generate.data3d() # list giving the locations to evaluate the basis functions x = dat$x # construct a set of basic B-splines for each dimension splines = default.splines(x) # construct assembled splines from splines list a = assemble(splines, x) # imagine there are 4 data obsevations we want to smooth # but that they can't be loaded into memory Ytilde = vector("list", 4) sum_ysq = numeric(4) # prepare each data set manually # notice the use of the assembled arguments so that # the splines are not "assembled" again for each data set for(i in seq_along(Ytilde)) { data = generate.data3d()$data3d Ytilde[[i]] = prepare(data, x = x, splines = splines, assembled = a) sum_ysq[i] = sum(data^2) } n = dim(data) p = create.prepared_list(assembled = a, x = x, Ytilde = Ytilde, sum_ysq = sum_ysq, n = n)# generate and prepare 3d data set.seed(9) dat = generate.data3d() # list giving the locations to evaluate the basis functions x = dat$x # construct a set of basic B-splines for each dimension splines = default.splines(x) # construct assembled splines from splines list a = assemble(splines, x) # imagine there are 4 data obsevations we want to smooth # but that they can't be loaded into memory Ytilde = vector("list", 4) sum_ysq = numeric(4) # prepare each data set manually # notice the use of the assembled arguments so that # the splines are not "assembled" again for each data set for(i in seq_along(Ytilde)) { data = generate.data3d()$data3d Ytilde[[i]] = prepare(data, x = x, splines = splines, assembled = a) sum_ysq[i] = sum(data^2) } n = dim(data) p = create.prepared_list(assembled = a, x = x, Ytilde = Ytilde, sum_ysq = sum_ysq, n = n)
evalargs
Create a default evalargs object based on
data. This is just a list of sequences.
If ni = dim(data)[i], then the sequence for
dimension i is seq(0, 1, len = ni).
default.evalargs(data)default.evalargs(data)
data |
A matrix or array-like object |
A list of equidistance sequences between 0 and 1
Joshua French
a = array(rnorm(10 * 11 * 12), dim = 10:12) default.evalargs(a)a = array(rnorm(10 * 11 * 12), dim = 10:12) default.evalargs(a)
Construct a list of hero_bsplines using the
default values suggested by Ruppert, Wand, and Carroll
(2003). Specifically, if
r = range(evalargs[[i]]) and
l = length(evalargs[[i]]), then Ruppert, Wand,
and Carroll (2003) suggest
nknots = min(ceiling(l/4), 35) and the function
returns the hero_bspline for that dimension as
bspline(r, nknots = nknots).
default.splines(evalargs)default.splines(evalargs)
evalargs |
A list of equidistant sequences. |
A list of hero_bsplines.
Joshua French
Ruppert, D., Wand, M. P., & Carroll, R. J. (2003). Semiparametric Regression. Cambridge University Press. <doi:10.1017/CBO9780511755453>
s1 = seq(0, 1, len = 10) s2 = seq(0, 1, len = 20) default.splines(list(s1, s2))s1 = seq(0, 1, len = 10) s2 = seq(0, 1, len = 20) default.splines(list(s1, s2))
P-spline difference penalty
diffpen(x, m = 2, sparse = TRUE)diffpen(x, m = 2, sparse = TRUE)
x |
A |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
A matrix or sparseMatrix-class object.
Joshua French
b = bspline(nbasis = 10) diffpen(b)b = bspline(nbasis = 10) diffpen(b)
enhance enhances the sandwich smoother by choosing
the optimal penalty value that minimizes the GCV
statistic. The optimx function
is used to do the optimization.
enhance( obj, par = rep(0, length(obj$n)), lower = rep(-20, length(par)), upper = rep(20, length(par)), method = "L-BFGS-B", control = list(), prepare = TRUE, loggcv = FALSE, ... )enhance( obj, par = rep(0, length(obj$n)), lower = rep(-20, length(par)), upper = rep(20, length(par)), method = "L-BFGS-B", control = list(), prepare = TRUE, loggcv = FALSE, ... )
obj |
A |
par |
a vector of initial values for the parameters for which optimal values are to be found. Names on the elements of this vector are preserved and used in the results data frame. |
lower, upper
|
Bounds on the variables for methods such as |
method |
The method to be used for optimization. The
default is |
control |
A list of control parameters. See ‘Details’. |
prepare |
A logical value. The default is |
loggcv |
A logical value indicating whether the log
of the GCV statistic should be used. Useful for very large
data sets. Default is |
... |
Additional arguments to pass to to the
|
Internally, the loglambda2gcv is used as
the objective function for the
optimx function. Many different
optimization methods are available. The default is
L-BFGS-B, which allows for constraints on the
parameters to optimize. Another excellent choice is the
nlminb algorithm, which also allows for parameter
constraints.
By default, a prepared_data object with
the optimal loglambda values that minimize the
GCV, along with an additional component,
results, that contains the optimization results.
Joshua French
# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) enhance(obj)# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) enhance(obj)
enhance.grid enhances the sandwich smoother by
choosing a optimal penalty value to lower the GCV
statistic. A grid search algorithm is utilized based on
the each row of par. The penalty values (assumed
to be on the log scale) are passed to the
loglambda2gcv function. If prepare
is TRUE, then obj is returned with the
penalty values that minimize the GCV statistic during the
grid search. Otherwise, the complete results of the grid
search are returned.
enhance.grid(obj, par, prepare = TRUE, loggcv = FALSE, ..., cl = NULL)enhance.grid(obj, par, prepare = TRUE, loggcv = FALSE, ..., cl = NULL)
obj |
A |
par |
A matrix-like object (i.e.,
|
prepare |
A logical value. The default is |
loggcv |
A logical value indicating whether the log
of the GCV statistic should be used. Useful for very large
data sets. Default is |
... |
Additional arguments to pass to to the
|
cl |
A cluster object created by |
By default, a prepared_* object with the
optimal loglambda values that minimize the GCV,
along with an additional component, results,
that contains the optimization results. Otherwise, the
complete results of the grid search.
Joshua French
# create b-splines b1 = bspline(nbasis = 10) b2 = bspline(nbasis = 12) # observed data locations x1 = seq(0, 1, len = 60) x2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(x1[i]-.5)^3)*cos(4*pi*x2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(x1, x2), list(b1, b2)) enhance.grid(obj, prepare = FALSE)# create b-splines b1 = bspline(nbasis = 10) b2 = bspline(nbasis = 12) # observed data locations x1 = seq(0, 1, len = 60) x2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(x1[i]-.5)^3)*cos(4*pi*x2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(x1, x2), list(b1, b2)) enhance.grid(obj, prepare = FALSE)
Enlarge the spatial domain of a
SpatialPolygons-class object. If
width isn't specified, then 10% of the maximum
distance between the points specified by the bounding
box is used. The
st_buffer function is used to enlarge
x.
enlarge(x, width, ...)enlarge(x, width, ...)
x |
A |
width |
The width to enlarge the study area. Distance from original geometry to include in the new geometry. Negative values are allowed. |
... |
Additional arguments to pass to
|
An object of class hero_enlarge. This is
simply a list with eborder (the enlarged
border), border (the border of the original
coordinates), and the width of the enlargement.
eborder and border are
SpatialPolygons-class objects.
Joshua French
# enlarge regular grid # create x and y defining square border x = seq(min(lon), max(lon), length = 60) y = seq(min(lat), max(lat), length = 80) border = border.grid(x, y) e = enlarge(border) plot(e) # create x and y defininging an irregular grid border2 = border.grid(lon, lat) e2 = enlarge(border2) plot(e2)# enlarge regular grid # create x and y defining square border x = seq(min(lon), max(lon), length = 60) y = seq(min(lat), max(lat), length = 80) border = border.grid(x, y) e = enlarge(border) plot(e) # create x and y defininging an irregular grid border2 = border.grid(lon, lat) e2 = enlarge(border2) plot(e2)
Generate two-dimensional data related to the f1 function
of Lu et al. (2012) (code from author). Define n =
c(60, 80). Then x[[i]] = (1:n[i])/n[i] -
1/2/n[i]. These are the observed data locations. For
i and j spanning the full length of each
element of x, mu2d[i, j] = sin(2 * pi *
(x[[1]][i] - .5) ^ 3) * cos(4 * pi * x[[2]][j]). Lastly,
data2d = mu2d + rnorm(prod(n)).
generate.data2d() generate_data2d() generateData2d() GenerateData2d()generate.data2d() generate_data2d() generateData2d() GenerateData2d()
A list with components x, mu2d, and
data2d. x is a list of sequences with
length 60 and 80. mu2d and data2d are
matrices of size 60 by 80.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
dat = generate.data2d()dat = generate.data2d()
Generate data related to Section 7.2 of Lu et al. (2012)
(code from author). Define n = c(128,
128, 24). Then x[[i]] = (1:n[i])/n[i] -
1/2/n[i]. These are the observed data locations. For
i, j, k spanning the full length
of each element of x, mu3d[i, j, k] =
x[[1]][i]^2 + x[[2]][j]^2 + x[[3]][k]^2. Lastly,
data3d = mu3d + 0.5 * rnorm(n[1] * n[2] * n[3]).
generate.data3d() generate_data3d() generateData3d() GenerateData3d()generate.data3d() generate_data3d() generateData3d() GenerateData3d()
A list with components x, mu3d, and
data3d. x is a list of sequences with
length 128, 128, and 24. mu3d and data3d
are arrays of size 128 by 128 by 24.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
dat = generate.data3d()dat = generate.data3d()
hero constructs a hero sandwich smoother based off
off a prepared data object coming from the
prepare function.
Subclasses are added (e.g., hero_numeric,
hero_matrix, hero_array, etc.) are added to
the returned object for plotting purposes.
A list is returned (and the data locations are not) for
hero.prepared_list. Each element of the list
contains the coefficients and fitted values (if fitted is
TRUE) for the respective data observation.
hero(x, ...) ## S3 method for class 'prepared_array' hero(x, ...) ## S3 method for class 'prepared_list' hero(x, ..., fitted = FALSE) ## S3 method for class 'prepared_matrix' hero(x, ...) ## S3 method for class 'prepared_numeric' hero(x, ...) ## S3 method for class 'prepared_sequential' hero( x, ..., export_list, export_fun = base::saveRDS, package = "base", call_args = list() ) ## S3 method for class 'prepared_starray' hero(x, ...) ## S3 method for class 'prepared_sts' hero(x, ...)hero(x, ...) ## S3 method for class 'prepared_array' hero(x, ...) ## S3 method for class 'prepared_list' hero(x, ..., fitted = FALSE) ## S3 method for class 'prepared_matrix' hero(x, ...) ## S3 method for class 'prepared_numeric' hero(x, ...) ## S3 method for class 'prepared_sequential' hero( x, ..., export_list, export_fun = base::saveRDS, package = "base", call_args = list() ) ## S3 method for class 'prepared_starray' hero(x, ...) ## S3 method for class 'prepared_sts' hero(x, ...)
x |
Data prepared via the |
... |
Mostly not implemented. |
fitted |
A logical value indicating whether the fitted values should be
computed. The default is |
export_list |
A vector or list whose elements tell |
export_fun |
A function that will write the results for each observation
to file using the names in |
package |
A character string indicating the approach to use for the
computations. The choices are |
call_args |
A named list providing relevant arguments to
|
A hero object with the smoothed data
(fitted), the estimated coefficients for the
basis functions (coefficients), and the
locations of the original data (x).
Joshua French.
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
French, Joshua P., and Piotr S. Kokoszka. "A sandwich smoother for spatio-temporal functional data." Spatial Statistics 42 (2021): 100413.
# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) obj = enhance(obj) sandmod = hero(obj) plot(sandmod)# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) obj = enhance(obj) sandmod = hero(obj) plot(sandmod)
See Details of bspline for additional
information about arguments.
knot.design( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE ) knot_design( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE ) knotDesign( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE ) KnotDesign( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE )knot.design( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE ) knot_design( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE ) knotDesign( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE ) KnotDesign( rangeval = 0:1, nbasis, nknots, norder = 4, extend = FALSE, interior = FALSE )
rangeval |
A numeric vector of length 2 defining the
interval over which the functional data object can be
evaulated. The default value is |
nbasis |
An integer specifying the number of
basis functions to construct. This is closely linked to
the number of knots ( |
nknots |
The number of *interior* knots. See Details. |
norder |
An integer specifying the order of the B-splines, which is one higher than their degree. The default is 4, which corresponds to cubic splines. |
extend |
Should the knots stop at the endpoints specified by |
interior |
A logical value specifying whether only interior knots should be returned. Default is |
An ascending sequence of univarite knot locations.
if (requireNamespace("fda", quietly = TRUE)) { b = fda::create.bspline.basis(nbasis = 10) # interior knots only bknots = b$params # should match knots = knot.design(nbasis = 10, interior = TRUE) all.equal(bknots, knots) }if (requireNamespace("fda", quietly = TRUE)) { b = fda::create.bspline.basis(nbasis = 10) # interior knots only bknots = b$params # should match knots = knot.design(nbasis = 10, interior = TRUE) all.equal(bknots, knots) }
A sequence of kronecker products
kronecker.seq(X, FUN = "*", make.dimnames = FALSE, ...) kronecker_seq(X, FUN = "*", make.dimnames = FALSE, ...) kroneckerSeq(X, FUN = "*", make.dimnames = FALSE, ...) KroneckerSeq(X, FUN = "*", make.dimnames = FALSE, ...)kronecker.seq(X, FUN = "*", make.dimnames = FALSE, ...) kronecker_seq(X, FUN = "*", make.dimnames = FALSE, ...) kroneckerSeq(X, FUN = "*", make.dimnames = FALSE, ...) KroneckerSeq(X, FUN = "*", make.dimnames = FALSE, ...)
X |
A list of numeric matrices or arrays |
FUN |
a function; it may be a quoted string. |
make.dimnames |
Provide dimnames that are the product of the
dimnames of |
... |
optional arguments to be passed to |
A matrix or array
x1 = matrix(rnorm(16), nrow = 4) x2 = matrix(rnorm(25), nrow = 5) x3 = matrix(rnorm(36), nrow = 6) x4 = matrix(rnorm(49), nrow = 7) p1 = x1 %x% x2 %x% x3 %x% x4 p2 = kronecker.seq(list(x1, x2, x3, x4)) all.equal(p1, p2)x1 = matrix(rnorm(16), nrow = 4) x2 = matrix(rnorm(25), nrow = 5) x3 = matrix(rnorm(36), nrow = 6) x4 = matrix(rnorm(49), nrow = 7) p1 = x1 %x% x2 %x% x3 %x% x4 p2 = kronecker.seq(list(x1, x2, x3, x4)) all.equal(p1, p2)
loglambda2gcv uses a vector of penalty values
to evaluate the GCV statistic for a
prepared_response object.
loglambda2gcv(loglambda, obj, loggcv = FALSE)loglambda2gcv(loglambda, obj, loggcv = FALSE)
loglambda |
A vector of penalty values (assumed to be on a natural logarithmic scale) for computing the GCV. |
obj |
A |
loggcv |
A logical value indicating whether the log of the
GCV statistic should be returned.
The default is |
Though this function can be used by the user, it is
basically an internal function used to find the
value of loglambda minimizing the GCV statistic.
The scalar GCV statistic
n1 = 10 b1 = bspline(nbasis = 10) x1 = seq(0, 1, len = n1) n2 = 20 x2 = seq(0, 1, len = n2) b2 = bspline(nbasis = 12) # construct "true" data mu = matrix(0, nrow = n1, ncol = n2) for(i in seq_len(n1)) { for(j in seq_len(n2)) { mu[i, j] = sin(2*pi*(x1[i]-.5)^3)*cos(4*pi*x2[j]) } } image(mu) # construct noisy data data = mu + rnorm(n1 * n2) x = list(x1, x2) splines = list(b1, b2) obj = prepare(data, x, splines) loglambda2gcv(c(0, 0), obj)n1 = 10 b1 = bspline(nbasis = 10) x1 = seq(0, 1, len = n1) n2 = 20 x2 = seq(0, 1, len = n2) b2 = bspline(nbasis = 12) # construct "true" data mu = matrix(0, nrow = n1, ncol = n2) for(i in seq_len(n1)) { for(j in seq_len(n2)) { mu[i, j] = sin(2*pi*(x1[i]-.5)^3)*cos(4*pi*x2[j]) } } image(mu) # construct noisy data data = mu + rnorm(n1 * n2) x = list(x1, x2) splines = list(b1, b2) obj = prepare(data, x, splines) loglambda2gcv(c(0, 0), obj)
Data related to the f1 function in
Lu et al. (2012). Define n1 = 60 and
x = seq_len(n1)/n1 - 1/2/n1. Similarly,
define n2 = 80 and
z = seq_len(n2)/n2 - 1/2/n2. The f1 function is defined as
sin(2 * pi * (x[i] - .5) ^ 3) * cos(4 * pi * z[j]),
where i in seq_along(x) and
j in seq_along(z). The result of this function
is stored in lutruef1. Using set.seed(3)
and adding rnorm(60 * 80) to lutruef1
results in lunoisyf1.
data(ludata)data(ludata)
The sequences x and z, the
lutruef1 data matrix with 60 rows and 80 columns,
and the lunoisyf1 data matrix with 60 rows and 80
columns.
hero_adjacent objectPlot a hero_adjacent object. x$nbrs is
used to construct a
sparseMatrix-class object. The
default behavior is to plot the sparse matrix using the
image function. However, if the
igraph package is installed, a graph is made using
graph_from_adjacency_matrix and
then plotted using plot.igraph.
## S3 method for class 'hero_adjacent' plot(x, ...)## S3 method for class 'hero_adjacent' plot(x, ...)
x |
A |
... |
Additional arguments passed to
|
coords = expand.grid(1:4, 1:4) a = adjacent(coords, digits = 1) plot(a)coords = expand.grid(1:4, 1:4) a = adjacent(coords, digits = 1) plot(a)
hero_bspline objectPlots basis functions specified by results of bspline.
## S3 method for class 'hero_bspline' plot(x, nderiv = 0L, type = "l", kcol = NULL, ...)## S3 method for class 'hero_bspline' plot(x, nderiv = 0L, type = "l", kcol = NULL, ...)
x |
An object of class |
nderiv |
An integer value specifying the derivative order of the B-splines. The default is 0. |
type |
character string (length 1 vector) or vector of 1-character
strings indicating the type of plot for each
column of |
kcol |
Color for vertical lines drawn at interior knots. Default is |
... |
Additional graphical parameters passed to |
x = bspline(nbasis = 10, extend = FALSE) plot(x) plot(x, nderiv = 1) plot(x, kcol = "grey") # plot vertical lines at knots # extend knots passed rangeval x2 = bspline(nbasis = 10, extend = TRUE) plot(x2, kcol = "grey") # compare to plot.fd if (requireNamespace("fda", quietly = TRUE)) { x3 = fda::create.bspline.basis(nbasis = 10) par(mfrow = c(2, 1)) plot(x, kcol = "grey") title("plot.hero_bspline") } plot(x3) title("plot.fd")x = bspline(nbasis = 10, extend = FALSE) plot(x) plot(x, nderiv = 1) plot(x, kcol = "grey") # plot vertical lines at knots # extend knots passed rangeval x2 = bspline(nbasis = 10, extend = TRUE) plot(x2, kcol = "grey") # compare to plot.fd if (requireNamespace("fda", quietly = TRUE)) { x3 = fda::create.bspline.basis(nbasis = 10) par(mfrow = c(2, 1)) plot(x, kcol = "grey") title("plot.hero_bspline") } plot(x3) title("plot.fd")
hero_enlarge objectPlot the enlarged and original border defined be a set of coordinates.
## S3 method for class 'hero_enlarge' plot(x, ..., blist = list(col = "grey"))## S3 method for class 'hero_enlarge' plot(x, ..., blist = list(col = "grey"))
x |
An object of class |
... |
Additional graphical parameters passed to the
plotting method for |
blist |
A list of additional graphical parameters
passed to the plotting method for |
b = border.grid(lon, lat) e = enlarge(b) plot(e)b = border.grid(lon, lat) e = enlarge(b) plot(e)
hero objectPlot the smoothed data produced by the
hero
function. The behavior of the function changes depending
on the subclass of the hero object. See Details.
## S3 method for class 'hero_matrix' plot(x, xlab = "", ylab = "", ...) ## S3 method for class 'hero_numeric' plot(x, xlab = "", ylab = "", type = "l", ...)## S3 method for class 'hero_matrix' plot(x, xlab = "", ylab = "", ...) ## S3 method for class 'hero_numeric' plot(x, xlab = "", ylab = "", type = "l", ...)
x |
An object of class |
xlab |
x-axis label |
ylab |
y-axis label |
... |
Additional graphical parameters passed to the relevant plotting function. See Details. |
type |
The plot type (when |
If x has subclass hero_numeric, then
the traditional plot function
is used to plot the smoothed data, with type = "l".
If x has subclass hero_matrix, then
image is used to plot the
smoothed data, or if the autoimage package is installed,
autoimage is used to
plot the smoothed data.
# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) obj = enhance(obj) sandmod = hero(obj) plot(sandmod)# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) obj = enhance(obj) sandmod = hero(obj) plot(sandmod)
hero_radspline
Plot a hero_radspline to compare the knots to the
observed data locations.
## S3 method for class 'hero_radspline' plot( x, blist = list(col = "grey"), glist = list(col = seq_along(x$grid) + 1, pch = seq_along(x$grid)), ... )## S3 method for class 'hero_radspline' plot( x, blist = list(col = "grey"), glist = list(col = seq_along(x$grid) + 1, pch = seq_along(x$grid)), ... )
x |
A |
blist |
A list to pass the plot method associated
with |
glist |
A list to pass the plot method associated
with |
... |
Additional arguments to pass the plot method
associated with |
If the default plotting styles for x$grid are to
be changed, the user can either choose a single
color/style that is replicated for each element of
x$grid or supply a vector which has length
matching length{x$grid}. See Examples.
Joshua French
border = border.grid(lon, lat) r = radspline(nknots = c(36, 36 * 4), border = border) # default color scheme plot(r) # change color and point styles of points, # and background of original domain plot(r, blist = list(col = "yellow"), glist = list(col = c("blue", "orange"), pch = 3:4))border = border.grid(lon, lat) r = radspline(nknots = c(36, 36 * 4), border = border) # default color scheme plot(r) # change color and point styles of points, # and background of original domain plot(r, blist = list(col = "yellow"), glist = list(col = c("blue", "orange"), pch = 3:4))
SpatialPolygons objectThis function takes a simple polygon and attempts to
convert it to a SpatialPolygons object.
This list is assumed to have components x and
y that define the boundary of the polygon.
poly2SpatialPolygons(x, ID = "border")poly2SpatialPolygons(x, ID = "border")
x |
A list with components |
ID |
The name of the resulting polygon. Default is
|
A SpatialPolygons object
Joshua French
angle = seq(0, 2 * pi, len = 100) poly = list(x = cos(angle), y = sin(angle)) plot(poly, type = "l", asp = 1) sppoly = poly2SpatialPolygons(poly) library(sp) plot(sppoly, axes = TRUE, asp = 1)angle = seq(0, 2 * pi, len = 100) poly = list(x = cos(angle), y = sin(angle)) plot(poly, type = "l", asp = 1) sppoly = poly2SpatialPolygons(poly) library(sp) plot(sppoly, axes = TRUE, asp = 1)
This function is an internal function to compute objects needed for fast implementation of the sandwich smoother. It is meant to be an internal function, so use this at your own risk.
precompute(B, P, m)precompute(B, P, m)
B |
A matrix of basis functions |
P |
A penalty matrix |
m |
Difference order of P-spline |
A list of needed objects
object = bspline(nbasis = 10) # sequence to evaluate evalarg = seq(0, 1, len = 11) # penalty matrix D = diffpen(object) P = Matrix::crossprod(D) B = predict(object, evalarg) stuff = precompute(B, P, m = 2)object = bspline(nbasis = 10) # sequence to evaluate evalarg = seq(0, 1, len = 11) # penalty matrix D = diffpen(object) P = Matrix::crossprod(D) B = predict(object, evalarg) stuff = precompute(B, P, m = 2)
hero objectPredict new values based on object produced by the
hero function.
## S3 method for class 'hero' predict(object, newB, ...)## S3 method for class 'hero' predict(object, newB, ...)
object |
A |
newB |
A vector or list containing the evaluated basis functions for the observations for which predictions are desired. |
... |
Not currently implemented. |
A matrix of the appropriate size
# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) obj = enhance(obj) sandmod = hero(obj) plot(sandmod) newb1 = predict(x1, newx = seq(0, 1, len = 100)) newb2 = predict(x2, newx = seq(0, 1, len = 100)) newB = list(newb1, newb2) p = predict(sandmod, newB = list(newb1, newb2))# create b-splines x1 = bspline(nbasis = 10) x2 = bspline(nbasis = 12) # observed data locations evalarg1 = seq(0, 1, len = 60) evalarg2 = seq(0, 1, len = 80) # construct "true" data mu = matrix(0, nrow = 60, ncol = 80) for(i in seq_len(60)) { for(j in seq_len(80)) { mu[i, j] = sin(2*pi*(evalarg1[i]-.5)^3)*cos(4*pi*evalarg2[j]) } } # construct noisy data data = mu + rnorm(60 * 80) obj = prepare(data, list(evalarg1, evalarg2), list(x1, x2)) obj = enhance(obj) sandmod = hero(obj) plot(sandmod) newb1 = predict(x1, newx = seq(0, 1, len = 100)) newb2 = predict(x2, newx = seq(0, 1, len = 100)) newB = list(newb1, newb2) p = predict(sandmod, newB = list(newb1, newb2))
hero_bspline objectPredicted values based on object created by
bspline.
## S3 method for class 'hero_bspline' predict(object, newx, nderiv = 0L, sparse = TRUE, ...)## S3 method for class 'hero_bspline' predict(object, newx, nderiv = 0L, sparse = TRUE, ...)
object |
A |
newx |
A numeric vector of values at which to evaluate the B-spline functions or derivatives. |
nderiv |
An integer value specifying the derivative order of the B-splines. The default is 0. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
... |
Not currently implemented. |
An matrix (or
Matrix-class object if
sparse = TRUE), where is the number of
values in newx and is the number of
basis functions in object. Each row gives the
predicted values of the basis functions for the
appropriate value of newx.
b = bspline(nbasis = 10) p = predict(b, newx = seq(0, 1, len = 101))b = bspline(nbasis = 10) p = predict(b, newx = seq(0, 1, len = 101))
hero_radspline
Predicted values based on object created by
radspline.
## S3 method for class 'hero_radspline' predict(object, newx, sparse = TRUE, longlat = FALSE, join = TRUE, ...)## S3 method for class 'hero_radspline' predict(object, newx, sparse = TRUE, longlat = FALSE, join = TRUE, ...)
object |
A |
newx |
A numeric matrix at which to evaluate the radial basis splines functions. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
longlat |
Use Euclidean ( |
join |
A logical value. |
... |
Not currently implemented. |
An matrix (or
Matrix-class object if
sparse = TRUE), where is the number of
rows in newx and is the number of
basis functions in object. Each row gives the
predicted values of each newx value evaluated
at each of the basis functions.
border = border.grid(lon, lat) r = radspline(nknots = c(36, 36 * 4), border = border) newx = cbind(c(lon), c(lat)) p = predict(r, newx)border = border.grid(lon, lat) r = radspline(nknots = c(36, 36 * 4), border = border) newx = cbind(c(lon), c(lat)) p = predict(r, newx)
A generic function to prepare various types of data. See the functions linked in See Also.
prepare(data, ...)prepare(data, ...)
data |
The data to prepare |
... |
Not implemented |
A prepared object
prepare.numeric,
prepare.matrix,
prepare.array,
prepare.sts,
prepare.starray
Sequentially prepare each observation for smoothing. It is assumed that each
observation resides in its own file and that do.call(import_fun,
list(import_list[i])) will import the data associated with observation
i into memory. The import_fun argument should be a function
after the style of readRDS, where the object can be
assigned a name once it is read in. The import_fun argument should NOT
be like load, where the object loaded has a preassigned
name.
prepare_sequential( import_list, import_fun = base::readRDS, x, splines, assembled, package = "base", call_args = list(), ... )prepare_sequential( import_list, import_fun = base::readRDS, x, splines, assembled, package = "base", call_args = list(), ... )
import_list |
A vector or list whose elements tell |
import_fun |
A function that will read each observation into memory
based on the elements of |
x |
The list of arguments at which to evaluate each
of the splines used to construct |
splines |
A list of spline-related objects. Each element of
|
assembled |
A list of |
package |
A character string indicating the package to use for the
computations. The choices are |
call_args |
A named list providing relevant arguments to the
|
... |
Not implemented |
A prepared_sequential object
Joshua P. French
prepare, mclapply,
pblapply, future_lapply,
mpi.applyLB
prepare.array prepares a data matrix for the
sandwich smooth. The dimensionality of data and
the length of x must match. Specifically,
length(dim(data)) must equal
length(x). The dimensionality of
data and the length of splines must match.
Specifically, length(dim(data)) must equal
length(splines).
## S3 method for class 'array' prepare(data, x, splines, m = 2, sparse = TRUE, ...)## S3 method for class 'array' prepare(data, x, splines, m = 2, sparse = TRUE, ...)
data |
A data array |
x |
A list of univariate, equidistant sequences. These should correspond to where the data are observed. Equidistant spacing between 0 and 1 is assumed if not supplied. See Details. |
splines |
A list of spline-related objects, e.g.,
produced by |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
... |
Not currently implemented. |
For a typical sandwich smooth, for data with
dimensions, Y[i1, i2, ...,id] is assumed to be
observed at position x[[1]][i1],
x[[2]][i2], ..., x[[d]][id].
Consequently, dim(data)[i] should equal
length(x[[i]]) for all i in
seq_len(d).
If x is not supplied, then
default.evalargs is used to create it
automatically.
If splines is not supplied, then a B-spline basis
is automatically created for each dimension using
default.splines.
A prepared_array object.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
bspline,
default.evalargs,
default.splines
# generate and prepare 3d data set.seed(9) dat = generate.data3d() obj = prepare(dat$data3d, x = dat$x)# generate and prepare 3d data set.seed(9) dat = generate.data3d() obj = prepare(dat$data3d, x = dat$x)
prepare.list prepares a list of data for the
sandwich smooth. The class of each element of
the list must be identical.
The dimensionality of data[[i]] and
the length of x must match. Specifically,
length(dim(data[[i]])) must equal
length(x). The dimensionality of
data[[i]] and the length of splines must match.
Specifically, length(dim(data[[i]])) must equal
length(splines). Note: If the splines
are preassembled, these can be passed using the argument
assembled so that this computation is not reperformed.
## S3 method for class 'list' prepare(data, x, splines, m = 2, sparse = TRUE, ...)## S3 method for class 'list' prepare(data, x, splines, m = 2, sparse = TRUE, ...)
data |
A list of |
x |
A list of values at which to evaluate the basis functions. See Examples and Details. |
splines |
A list of spline objects
( |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
... |
Not currently implemented. |
This function applies the functions
prepare.numeric,
prepare.matrix, and
prepare.array to each element of the list,
so relevant restrictions in the arguments may be
found there.
A prepared_list object.
Joshua French.
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
prepare.numeric, prepare.matrix,
prepare.array
# generate and prepare 3d data set.seed(9) dat = lapply(1:3, function (i) generate.data3d()) x = dat[[1]]$x data = lapply(dat, getElement, name = "data3d") obj = prepare(data, x = x) h = hero(obj)# generate and prepare 3d data set.seed(9) dat = lapply(1:3, function (i) generate.data3d()) x = dat[[1]]$x data = lapply(dat, getElement, name = "data3d") obj = prepare(data, x = x) h = hero(obj)
prepare.matrix prepares a data matrix for the
sandwich smooth. The dimensionality of data and
the length of x must match. Specifically,
length(dim(data)) must equal
length(x). The dimensionality of
data and the length of splines must match.
Specifically, length(dim(data)) must equal
length(splines).
## S3 method for class 'matrix' prepare( data, x, splines, m = 2, sparse = TRUE, spdiffpen = TRUE, digits = 1, sts = FALSE, ... )## S3 method for class 'matrix' prepare( data, x, splines, m = 2, sparse = TRUE, spdiffpen = TRUE, digits = 1, sts = FALSE, ... )
data |
A data matrix. |
x |
A list of values at which to evaluate the basis functions. See Examples and Details. |
splines |
A list of spline objects
( |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
spdiffpen |
A logical value indicating whether
|
digits |
The number of digits to use when applying
|
sts |
A logical value indicating whether |
... |
Not currently implemented. |
For a typical sandwich smooth (sts = FALSE),
for two-dimensional data, data[i, j] is assumed
to be observed at position x[[1]][i],
x[[2]][j]. If the data are a spatial time series,
then the first dimension is assumed to refer to space,
and the second dimension to time. In that case,
data[i, j] is assumed
to be observed at location x[[1]][i, ] and time
x[[2]][j].
If sts = TRUE, then x[[1]] should be a
matrix of spatial coordinates, with each row
corresponding to a location, and x[[2]] should
be a vector with the observation times.
If x is not supplied, then
default.evalargs is used to create it
automatically. This is only valid when
sts = FALSE.
If splines is not supplied, then a B-spline basis
is automatically created for each dimension using
default.splines. This is only valid when
sts = FALSE.
A prepared_matrix object.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
bspline, radspline,
diffpen, spdiffpen,
default.evalargs,
default.splines
# prepare Lu et al. (2012) noisy f1 data data(ludata) obj = prepare(lunoisyf1, x = list(x, z)) h = hero(obj) # precompute some stuff splines = default.splines(list(x, z)) l = assemble(splines, x = list(x, z)) obj2 = prepare(lunoisyf1, x = list(x, z), splines = splines, assembled = l) h2 = hero(obj2) all.equal(h, h2)# prepare Lu et al. (2012) noisy f1 data data(ludata) obj = prepare(lunoisyf1, x = list(x, z)) h = hero(obj) # precompute some stuff splines = default.splines(list(x, z)) l = assemble(splines, x = list(x, z)) obj2 = prepare(lunoisyf1, x = list(x, z), splines = splines, assembled = l) h2 = hero(obj2) all.equal(h, h2)
prepare.vector prepares a data vector for the
sandwich smooth. Unlike the other prepare.*
functions, x and splines do not
need to be lists since the data are 1-dimensional.
## S3 method for class 'numeric' prepare(data, x, splines, m = 2, sparse = TRUE, ...)## S3 method for class 'numeric' prepare(data, x, splines, m = 2, sparse = TRUE, ...)
data |
A numeric data vector |
x |
A sequence of equidistant values corresponding to where the data are observed. Equidistant spacing between 0 and 1 is assumed if not supplied. See Details. |
splines |
A spline-related object, e.g.,
produced by |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
... |
Not currently implemented. |
If x is not supplied and n is
the length(data), then the function automatically
sets x = seq(0, 1, length = n).
If splines is not supplied, and n is the
length(data), then the function automatically sets
splines = bspline(range(x), nknots = min(ceiling(n/4), 35)).
A prepared_numeric object.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
Ruppert, D., Wand, M. P., & Carroll, R. J. (2003). Semiparametric Regression. Cambridge University Press. <doi:10.1017/CBO9780511755453>
bspline,
default.evalargs,
default.splines
# create data n = 160 x = seq(0, 4 * pi, len = n) # "true" data mu = sin(x) # plot true data plot(x, mu, type = "l") # construct noisy data set.seed(4) data = mu + rnorm(n) # construct spline splines = bspline(c(0, 4 * pi), nknots = 20) # prepare/enhance data obj = prepare(data, x, splines) obj = enhance(obj) sandmod = hero(obj) plot(sandmod, ylim = range(data), lty = 2) lines(x, data, col = "lightgrey") lines(x, mu) legend("bottomleft", legend = c("smoothed", "true", "observed"), lty = c(2, 1, 1), col = c("black", "black", "grey"))# create data n = 160 x = seq(0, 4 * pi, len = n) # "true" data mu = sin(x) # plot true data plot(x, mu, type = "l") # construct noisy data set.seed(4) data = mu + rnorm(n) # construct spline splines = bspline(c(0, 4 * pi), nknots = 20) # prepare/enhance data obj = prepare(data, x, splines) obj = enhance(obj) sandmod = hero(obj) plot(sandmod, ylim = range(data), lty = 2) lines(x, data, col = "lightgrey") lines(x, mu) legend("bottomleft", legend = c("smoothed", "true", "observed"), lty = c(2, 1, 1), col = c("black", "black", "grey"))
starray for sandwich smoothprepare.starray prepares a spatio-temporal
array for the sandwich smooth.
## S3 method for class 'starray' prepare(data, x, y, times, rs, bs, m = 2, sparse = TRUE, spdiffpen = TRUE, ...)## S3 method for class 'starray' prepare(data, x, y, times, rs, bs, m = 2, sparse = TRUE, spdiffpen = TRUE, ...)
data |
An |
x |
A vector or matrix of x coordinates. See Details. |
y |
A vector or matrix of y coordinates. See Details. |
times |
The vector of times at which the data were observed. |
rs |
A |
bs |
A |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
spdiffpen |
A logical value indicating whether
|
... |
Not currently implemented. |
A prepared_starray object.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
# construct basis functions border = border.grid(lon, lat) rs = radspline(nknots = 36, poverlap = 3, border = border, longlat = TRUE) bs = bspline(c(1, 30), nbasis = 6) data = starray(tasmax) p = prepare(data, x = lon, y = lat, times = 1:30, rs = rs, bs = bs)# construct basis functions border = border.grid(lon, lat) rs = radspline(nknots = 36, poverlap = 3, border = border, longlat = TRUE) bs = bspline(c(1, 30), nbasis = 6) data = starray(tasmax) p = prepare(data, x = lon, y = lat, times = 1:30, rs = rs, bs = bs)
starray for sandwich smoothprepare.starray prepares a spatio-temporal
array for the sandwich smooth.
## S3 method for class 'sts' prepare( data, coords, times, rs, bs, m = 2, sparse = TRUE, spdiffpen = TRUE, ... )## S3 method for class 'sts' prepare( data, coords, times, rs, bs, m = 2, sparse = TRUE, spdiffpen = TRUE, ... )
data |
An |
coords |
A two-dimensional matrix-like object with non-NULL dimensions. |
times |
The vector of times at which the data were observed. |
rs |
A |
bs |
A |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
spdiffpen |
A logical value indicating whether
|
... |
Not currently implemented. |
A prepared_sts object.
Joshua French. Based off code by Luo Xiao (see References).
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
# construct basis functions border = border.grid(lon, lat) rs = radspline(nknots = 36, poverlap = 3, border = border, longlat = TRUE) bs = bspline(c(1, 30), nbasis = 6) splines = list(rs, bs) data = as.sts(tasmax) p = prepare(data, coords = cbind(c(lon), c(lat)), times = 1:30, rs = rs, bs = bs)# construct basis functions border = border.grid(lon, lat) rs = radspline(nknots = 36, poverlap = 3, border = border, longlat = TRUE) bs = bspline(c(1, 30), nbasis = 6) splines = list(rs, bs) data = as.sts(tasmax) p = prepare(data, coords = cbind(c(lon), c(lat)), times = 1:30, rs = rs, bs = bs)
radspline specifies a set of radial basis splines.
nknots is the approximate number of knots to
sample in the (usually) enlarged study area. If
eborder is not provided, then eborder is
automatically constructed by enlarging the border
object using the enlarge function and
width. See Details for additional information
about sampling the knot locations.
radspline( nknots, border, poverlap = 2, k = 2, width, type = "hexagonal", longlat = FALSE, eborder, ... )radspline( nknots, border, poverlap = 2, k = 2, width, type = "hexagonal", longlat = FALSE, eborder, ... )
nknots |
The approximate number of knots locations. Can be a vector of positive integers for successive samplings. See Details. |
border |
A |
poverlap |
The proportional amount of overlap (>=1) beyond the nearest neighbor knots. Default is 2. |
k |
The order of the Wendland covariance function. |
width |
The width for the border enlargement. |
type |
The sampling type for
|
longlat |
A logical value indicating whether Great
Circle distances should be used ( |
eborder |
A |
... |
Additional arguments passed to
|
The spsample function is used to
"automatically" select the knot locations within
eborder. nknots corresponds to the
n argument in that function. A hexagonal sampling
scheme is used by default, but other options are
available.
Great circle distance IS NOT used in sampling from the
regular grid. This is computationally expensive, so it
has not been implemented. Great circle distance is only
used when the constructed hero_radspline is
evaluated (and longlat = TRUE).
A hero_radspline object.
border = border.grid(lon, lat) r = radspline(nknots = c(36, 36 * 4), border = border) # default color scheme plot(r) # change color and point styles of points, # and background of original domain plot(r, blist = list(col = "yellow"), glist = list(col = c("blue", "orange"), pch = 3:4))border = border.grid(lon, lat) r = radspline(nknots = c(36, 36 * 4), border = border) # default color scheme plot(r) # change color and point styles of points, # and background of original domain plot(r, blist = list(col = "yellow"), glist = list(col = c("blue", "orange"), pch = 3:4))
A rotation of the H-transform of the array a by a
matrix x.
rh(x, a, transpose = FALSE)rh(x, a, transpose = FALSE)
x |
A matrix-like object. See Details. |
a |
An d-dimensional array |
transpose |
A logical value. The Default is
|
x should be matrix-like. This function has been
tested when x is a matrix object or a
Matrix.
Assuming a is of size , then x is of size .
A rotated, h-transformed array
Joshua French. Based off code by Luo Xiao (see References).
Currie, I. D., Durban, M. and Eilers, P. H. (2006), Generalized linear array models with applications to multidimensional smoothing. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68: 259-280. <doi:10.1111/j.1467-9868.2006.00543.x>
Xiao, L. , Li, Y. and Ruppert, D. (2013), Fast bivariate P-splines: the sandwich smoother. J. R. Stat. Soc. B, 75: 577-599. <doi:10.1111/rssb.12007>
dim = c(10:12) # construct random array a = array(rnorm(prod(dim)), dim = dim) # construct random matrix x = matrix(rnorm(15 * dim[1]), nrow = 15) rhxa = rh(x, a)dim = c(10:12) # construct random array a = array(rnorm(prod(dim)), dim = dim) # construct random matrix x = matrix(rnorm(15 * dim[1]), nrow = 15) rhxa = rh(x, a)
rh sequentiallyrh.seq sequentially applies the rh
function to a. Specifically, if the length of
x is d, then rh.seq(x, a) is
equivalent to rh(x[[d]], rh(x[[d - 1]], ..., rh(x[[2]], rh(x[[1]], a))..)).
rh.seq(x, a, transpose = FALSE) rh_seq(x, a, transpose = FALSE) rhSeq(x, a, transpose = FALSE) RhSeq(x, a, transpose = FALSE)rh.seq(x, a, transpose = FALSE) rh_seq(x, a, transpose = FALSE) rhSeq(x, a, transpose = FALSE) RhSeq(x, a, transpose = FALSE)
x |
A list of matrix-like objects |
a |
A matrix-like object (with dimensions) |
transpose |
A logical value. The Default is
|
A matrix or Matrix-class.
# generate x, a x = list(matrix(rnorm(100), nrow = 10), matrix(rnorm(100), nrow = 10)) a = matrix(rnorm(100), nrow = 10) # three equivalent forms rhs1 = rh.seq(x, a) rhs2 = rh(x[[2]], rh(x[[1]], a)) rhs3 = x[[1]] %*% a %*% t(x[[2]]) # check equality all.equal(rhs1, rhs2) all.equal(rhs1, rhs3)# generate x, a x = list(matrix(rnorm(100), nrow = 10), matrix(rnorm(100), nrow = 10)) a = matrix(rnorm(100), nrow = 10) # three equivalent forms rhs1 = rh.seq(x, a) rhs2 = rh(x[[2]], rh(x[[1]], a)) rhs3 = x[[1]] %*% a %*% t(x[[2]]) # check equality all.equal(rhs1, rhs2) all.equal(rhs1, rhs3)
spdiffpen computes the mth order spatial
difference penalty for a set of coordinates.
spdiffpen(coords, m = 1, sparse = TRUE, longlat = FALSE, digits = 1)spdiffpen(coords, m = 1, sparse = TRUE, longlat = FALSE, digits = 1)
coords |
A two-dimensional matrix-like object with non-NULL dimensions. |
m |
A positive integer indicating order of the difference penalty. |
sparse |
A logical value indicating if the result
should be a sparse version of the
|
longlat |
A logical value indicating whether Great
Circle distances should be used ( |
digits |
The number of digits to use when applying
|
adjacent is used to determine the
first-order neighbors of each point in coords. The
difference penalties are then successively determined
from that.
If sparse = TRUE, a
sparseMatrix-class Matrix is
returned when the penalty matrix is relatively sparse
(typically, at least half the entries are zero).
Otherwise, something of the more general
Matrix-class is returned.
A matrix or sparseMatrix-class object.
coords = expand.grid(1:4, 1:4) # first order difference penalty d1 = spdiffpen(coords, digits = 1) # second order difference penalty d2 = spdiffpen(coords, m = 2, digits = 1) # third order difference penalty d3 = spdiffpen(coords, m = 3, digits = 1)coords = expand.grid(1:4, 1:4) # first order difference penalty d1 = spdiffpen(coords, digits = 1) # second order difference penalty d2 = spdiffpen(coords, m = 2, digits = 1) # third order difference penalty d3 = spdiffpen(coords, m = 3, digits = 1)
The maximum daily surface air temperature (C) for the time period January 1, 1971 through January 30, 1971 for the ECP2-GFDL computer generated data made available through the North American Regional Climate Change Assessment Program (NARCCAP).
data(tasmax)data(tasmax)
Matrices lon and lat and array
tasmax.
Mearns, L.O., et al., 2007, updated 2014. The North American Regional Climate Change Assessment Program dataset, National Center for Atmospheric Research Earth System Grid data portal, Boulder, CO. Data downloaded 2018-06-13. <doi:10.5065/D6RN35ST>.
Mearns, L. O., W. J. Gutowski, R. Jones, L.-Y. Leung, S. McGinnis, A. M. B. Nunes, and Y. Qian. "A regional climate change assessment program for North America." EOS, Vol. 90, No. 36, 8 September 2009, pp. 311-312. <doi:10.1029/2009EO360002>.
The maximum daily surface air temperature (C) of land locations for the time period January 1, 2041 through January 30, 1941 for the WRFG-CGCM3 computer generated data made available through the North American Regional Climate Change Assessment Program (NARCCAP). The non-land locations are specified as NA.
data(wrfg_cgcm3_tasmax)data(wrfg_cgcm3_tasmax)
Matrices wrfg_lon and wrfg_lat and array
wrfg_cgcm3_tasmax.
Mearns, L.O., et al., 2007, updated 2014. The North American Regional Climate Change Assessment Program dataset, National Center for Atmospheric Research Earth System Grid data portal, Boulder, CO. Data downloaded 2018-06-13. <doi:10.5065/D6RN35ST>.
Mearns, L. O., W. J. Gutowski, R. Jones, L.-Y. Leung, S. McGinnis, A. M. B. Nunes, and Y. Qian. "A regional climate change assessment program for North America." EOS, Vol. 90, No. 36, 8 September 2009, pp. 311-312. <doi:10.1029/2009EO360002>.