fields (version 1.7.2)

image.cov: Exponential,Gaussian, "power" and Matern covariance functions for 2-d gridded locations

Description

Given two sets of locations defined on a 2-d grid efficiently multiplies a cross covariance with a vector.

Usage

exp.image.cov(ind1, ind2, Y, cov.obj = NULL, setup = FALSE, grid, ...) 
matern.image.cov(ind1, ind2, Y, cov.obj = NULL, setup = FALSE, grid,
M=NULL,N=NULL,...)

Value

A vector that is the multiplication of the cross covariance matrix with the vector Y.

Details

This function was provided to do fast computations for large numbers of spatial locations and supports the conjugate gradient solution in krig.image. In doing so the observations can be irregular spaced but their coordinates must be 2-dimensional and be restricted to grid points. (The function as.image will take irregular, continuous coordinates and overlay a grid on them.)

Returned value: If ind1 and ind2 are matrices where nrow(ind1)=m and nrow(ind2)=n then the cross covariance matrix, Sigma is an mXn matrix (i,j) element is the covariance between the grid locations indexed at ind1[i,] and ind2[j,]. The returned result is Sigma%*%Y. Note that one can always recover the coordinates themselves by evaluating the grid list at the indices. e.g. cbind( grid$x[ ind1[,1]], grid$y[ind1[,2])) will give the coordinates associated with ind1. Clearly it is better just to work with ind1! Functional Form: Following the same form as exp.cov and matern.cov for irregular locations, the covariance is defined as phi( D.ij) where D.ij is the Euclidean distance between x1[i,] and x2[j,] but having first been scaled by theta. Specifically,

D.ij = sqrt( sum.k (( x1[i,k] - x2[j,k]) /theta[k])**2 ).

For the exponential phi(u) = exp( -u**p) and for the Matern phi(u)= matern( u, smoothness). (See the help file on matern.cov for details.) Note that if theta is a scalar then this defines an isotropic covariance function.

Implementation: This function does the multiplication on the full grid efficiently by a 2-d FFT. The irregular pattern in Y is handled by padding with zeroes and once that multiplication is done only the appropriate subset is returned.

As an example assume that the grid is 100X100 let big.Sigma denote the big covariance matrix among all grid points ( If the parent grid is 100x100 then big.Sigma is 10K by 10K !) Here are the computing steps:

temp<- matrix( 0, 100,100)

temp[ ind2] <- Y

temp2<- big.Sigma%*% temp

temp2[ind1]

Notice how much we pad with zeroes or at the end throw away! Here the matrix multiplication is effected through convolution/FFT tricks to avoid creating and multiplying big.Sigma explicitly. It is often faster to multiply the regular grid and throw away the parts we do not need then to deal directly with the irregular set of locations.

Confusion: In this entire discussion Y is treated as vector. However is one has complete data then Y can also be interpreted as a image matrix conformed to correspond to spatial locations. See the last example for this distinction.

See Also

exp.cov, matern.cov,smooth.2d, as.image, krig.image

Examples

# multiply 2-d isotropic exponential with theta=4 by a random vector 
junk<- matrix(rnorm(100*100), 100,100)
cov.obj<- exp.image.cov( setup=TRUE, grid=list(x=1:100,y=1:100), theta=8)
result<-  exp.image.cov(Y=junk,cov.obj=cov.obj)

image( matrix( result, 100,100))
# do it again, no setup needed 
#  junk2<- matrix(rnorm(100**2, 100,10))
#  result2<-  exp.image.cov(Y=junk2, cov.obj=cov.obj)



# generate a grid and set of indices based on discretizing the locations
# in the precip dataset
data(RMprecip)
out<-as.image( RMprecip$y, x= RMprecip$x)
ind1<- out$ind
grid<- list( x= out$x, y=out$y)
# there are 750 gridded locations 

# setup to create cov.obj for exponential covariance with range= 1.25
# p=1 is the default
cov.obj<- exp.image.cov( setup=TRUE, grid=grid, theta=1.25) 

# multiply covariance matrix by an arbitrary vector
junk<-  rnorm(750)
result<- exp.image.cov( ind1, ind1, Y= junk,cov.obj=cov.obj)

# The brute force way would be  
# result<- exp.cov( RMprecip$x, RMprecip$x,theta=1.25, C=junk)


# evaluate the covariance between all grid points and the center grid point
Y<- matrix(0,cov.obj$m, cov.obj$n)
Y[32,32]<- 1
result<- exp.image.cov( Y= Y,cov.obj=cov.obj)
# covariance surface with respect to the grid point at (32,32)
# 
# reshape "vector" as an image
temp<-  matrix( result, cov.obj$m,cov.obj$n)
image.plot(cov.obj$grid$x,cov.obj$grid$y, temp)
# or persp( cov.obj$grid$x,cov.obj$grid$y, temp) 

# check out the Matern. 
cov.obj<- matern.image.cov( setup=TRUE, grid=grid, theta=1.25, 
smoothness=2)
Y<- matrix(0,64,64)
Y[16,16]<- 1
result<- matern.image.cov( Y= Y,cov.obj=cov.obj)
temp<-  matrix( result, cov.obj$m,cov.obj$n)
image.plot( cov.obj$grid$x,cov.obj$grid$y, temp)
# Note we have centered at the location (16,16) for this case