Radial.Basis: Two dimensional radial and tensor basis functions

In LatticeKrig: Multi-Resolution Kriging Based on Markov Random Fields

Two dimensional radial and tensor basis functions

Description

Two dimensional radial basis and tensor functions that are used as the building blocks for the LatticeKrig multiresolution sparse basis matrices.

Usage

Radial.basis(x1, centers, basis.delta, max.points = NULL,
                 mean.neighbor = 50, BasisFunction =
                 "WendlandFunction", distance.type = "Euclidean",
                 verbose = FALSE)

Tensor.basis(x1, centers, basis.delta, max.points = NULL, mean.neighbor = 50, 
   BasisFunction = "WendlandFunction", distance.type = "Euclidean",
   verbose = FALSE) 

WendlandFunction(d)    

triWeight(d)

CubicBSpline(d)

Arguments

x1	A matrix of locations to evaluate the basis functions. Each row of x1 is a location. For spherical distances these locations are assumed to be in degrees.
centers	A matrix specifying the basis function centers. For spherical distances these locations are assumed to be in degrees.
d	A vector of distances.
basis.delta	A vector of scale parameters for the basis functions.
max.points	Maximum number of nonzero entries expected for the returned matrix.
distance.type	A text string indicate type distance to use between spatial locations when evaluating the basis functions. Default is "Euclidean". Other choices when locations are in degrees ( e.g. longitude and latitude) are "Chordal" and "GreatCircle" with the default units being miles. See Details below how to change the radius that is used.
mean.neighbor	Average number of centers that are within delta of each x1 location. For centers on a regular grid this is often easy to estimate.
BasisFunction	A function that will take a non-negative argument and be zero outside [0,1]. This is applied to distance(s) to generate the basis functions. For tensor basis functions, the function is applied to the distance components for each dimension.
verbose	Print out debugging information if TRUE.

Details

Radial.basis and Tensor.basis These function finds the pairwise distances between the points x1 and centers and evaluates the function RadialBasisFunction at these distances scaled by delta. In most applications delta is constant, but a variable delta could be useful for lon/lat regular grids. The Wendland function is for 2 dimensions and smoothness order 2. See WendlandFunction for the polynomial form. This code has a very similar function to the fields function wendland.cov.

In pseudo R code for delta a scalar Radial.basis with the Wendland "shape" evaluates as

  BigD<- rdist( x1,centers)
  WendlandFunction(BigD/basis.delta)

and is comparable to

Radial.basis(x1, centers,basis.delta)

The actual code uses a FORTRAN subroutine to search over distances less than delta and also returns the matrix in sparse format.

Part of the flexibility in these function is to use different distance functions (metrics). See help on LKDistance for the S4 method to find distance and LKDist for the lower level functions that actually do the work. To change the radius used for the spherical distances one add a Radius attribute to the text string. For example to use kilometers for great circle distance ( R approximately 6371)

   dtype <- "GreatCircle" 
   attr( dtype, which="Radius") <- 6371 
   

dtype is still a text string but has the extra bit of information attached to it. Now use distance.type = dtype for the distance type argument.

The function Tensor.basis has similar function as the radial option. The main difference is that a slightly different distance function is used to return the component distances for each dimension. In pseudo R code for delta a scalar and for just two dimensions Tensor.basis evaluates as

  BigD1<- rdist( x1[,1],centers[,1])
  BigD2<- rdist( x1[,2],centers[,2])
  WendlandFunction(BigD1/basis.delta) *WendlandFunction(BigD1/basis.delta)

The higher level function LKrig.basis paired with an LKinfo object to define the centers returns a sparse matrix with the basis functions evaluated at a set of locations.

If normalize = FALSE in the LKinfo object then the basis functions are evaluated exactly as described above. If normalize = TRUE there is an additional set of weights that are computed so that the LatticeKrig model when viewed as a Gaussian process has constant marginal variances.

WendlandFunction, triWeight, CubicBSpline These functions have limited support and will be zero for d greater than one.

Wendland (\phi(d)) for d \le 1

\phi(d) = (1 - d)^6 (35d^2 + 18 d + 3))/3

and zero otherwise.

triWeight (\phi(d)) for d \le 1)

\phi(d) = (1 - d^2)^3

CubicBSpline (\phi(d)) for d \le 1)

A piecewise cubic polynomial over the four intervals

\{ [-1,-1/2],[-1/2,0],[0,-1/2],[1/2,1] \}

the individual cubic polynomials enforce zero boundary conditions for the function, first and second derivaitves at -1 and 1 and continuity of the function, first and second derivatives at the three interior knot points (\{ -1/2, 0, -1/2 } )

Value

For Radial.basis and Tensor.basis a matrix in sparse format with number of rows equal to nrow(x1) and columns equal to nrow(center) and value are the basis functions evaluated at these combinations For WendlandFunction and triweight a vector with the values for this function.

Author(s)

Doug Nychka

See Also

LKrig.basis, LKDist

Examples

set.seed(12)
x<- cbind( runif(100), runif(100))
center<- expand.grid( seq( 0,1,,5), seq(0,1,,5))
# coerce to matrix
center<- as.matrix(center)

  PHI1<- Radial.basis(x, center, basis.delta = .5)
  PHI2<- Tensor.basis( x, center, basis.delta = .5 )
# similarity of radial and tensor product forms  
  plot( c(0,1.1), c(0,1.1), type="p")
  for( k in 1:25){
	points( PHI1[,k], PHI2[,k])
	}
	
# LKrig with a different radial basis function. 
# 
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find locations that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]
  obj<- LKrig(x,y,NC=30,nlevel=1, alpha=1, lambda=.01, a.wght=5)
    
  obj1<- LKrig(x,y,NC=30,nlevel=1, alpha=1, 
    lambda=.01, a.wght=5, BasisFunction="triWeight", overlap=1.8)
    
# It is usually better to create the LKinfo object first and then call LKrig
LKinfo <- LKrigSetup( x,NC=30,nlevel=1, alpha=1, 
            lambda=.01, a.wght=5, BasisFunction="triWeight", overlap=1.8 )
obj1B<- LKrig(x,y, LKinfo=LKinfo)  
#######################################################
###  a radial basis on the sphere
#######################################################
# icosohedral grid at the second generation 
allLevels<- IcosahedronGrid(2)
# extract second level
centers3d<- allLevels[[2]]
centers<- toSphere(centers3d)
# take a look
plot( centers, xlim=c(-180,180), ylim=c( -90, 90), pch=16)
lonlat<- make.surface.grid( list( x= seq( -180,180,,80), y= seq( -90,90,,40)))
delta<- 5000 # default distance is in miles 
# evaluate all basis functions on the grid
bigX<- Radial.basis( lonlat, centers, basis.delta= delta, 
distance.type="GreatCircle")
# look at the 12th basis function 
b12<- bigX[,12]
image( as.surface( lonlat, b12), col=c( "white",tim.colors(256)), add=TRUE,
xlim="longitude", ylim="latitude")
points( centers,  pch=16, col="magenta", cex=2)
#######################################################
###  check of the cubicBspline
#######################################################
xgrid<- seq( -1.2,1.2,1e-3)
y<- CubicBSpline( xgrid)

xgrid<- seq( -1.2,1.2,.005)
y<- CubicBSpline( xgrid)
plot( xgrid, y, type="l")
rug( c( -1, -.5,0,.5,1),col="blue", lwd=4)

#If piecewise cubic third differences should be step functions
# with jumps at the knots
plot( xgrid[-(1:3)], diff( y,differences=3)/(xgrid[2]- xgrid[1])^3, type="l")
rug( c( -1, -.5,0,.5,1),col="blue", lwd=4)

LatticeKrig documentation built on May 30, 2026, 5:07 p.m.