Wimage.cov: Functions for W-transform based covariance models

Description

These functions are designed for 2-d Gaussian random fields on a large regular grid. They require a sparse or block diagonal model for the covariances among the wavelet coefficients. Wimage.cov multiplies a vector with the implied covariance matrix and Wimage.sim generates a random field with the implied covariance.

Usage

Wimage.cov(Y = NULL, ind = NULL, H.obj, find.row = FALSE)
Wimage.sim( H.obj))

Arguments

Value

An mXn matrix with the same extent as the grid defined by H.obj.

Details

Note that Wimage.cov only returns the product of a vector times the covariance matrix. This is usually all that is really needed and finesses the memeory problmes of dealing with very large covaraince matrices. But see the last example for a loop to extract a submatrix of the covariance.

The way 2-d wavelet basis functions are indexed can be confusing. See the help files for Wimage.info and Wtransform.image for more background. For these functions, ind can either be a single vector and the index refers to the grid points in column stacked (or ``vec'') format.

If ind has a two columns these refer to the row/column of the grid. Currently only one sparse representation for H is implemented, a block/diagonal strategy. The covariance has the form Wi*H *H * Wi.t. Where Wi is mnXmn the matrix of wavelet basis functions with each column being a basis evaluated on the gird points and in stacked column format. Wi.t is its transpose. H is represented as partitioned matrix where H0 is a full matrix and the remaining rows only have a diagonal term. If ind0 denotes the indices for the elements of H0 then in R notation:

H0 = H[ ind0, ind0]

To simplify the coding the diagonal elements are represented as an mXn matrix with those locations corresponding to H0 set to 1. Accordingly, in R notation

H1= diag(H), H1[ind0] <- 1

With this representation, multiplying H by a vector u becomes

u[ind0]<- H0%*% c( u[ind0])

u*H1

Note that u is assumed to be an mXn matrix based on the grid dimensions.

Under the assumptions the grid ( or image) is mXn the components for H.obj are m , n, cut.min, specifying the coarsest level of wavelet basis functions, H0, ind, and H1.

Examples

Run this code

# This example  creates a block/diagonal for H that is close to a 
# Matern covariance
# The example while lengthy is a practical run through of the process 
# of creating a reasonable sparse representation for H. 

M<- 16 # a 16X16 grid of locations
MM<- M**2
x<- 1:M
y<- 1:M
grid<- list( x=x, y=y)
cut.min<- 4
# H0 is father and mothers at first level
# with cut.min=4

# get indices associated with H0
I<-rep( 1:cut.min, cut.min)
J<- sort( I)
Wimage.s2i(I, J, flavor=0, level=0, cut.min=cut.min, m=M, n=M, mat=FALSE)-> ind0
Wimage.s2i(I, J, flavor=1, level=1, cut.min=cut.min, m=M, n=M, mat=FALSE)-> ind1
Wimage.s2i(I, J, flavor=2, level=1, cut.min=cut.min, m=M, n=M, mat=FALSE)-> ind2
Wimage.s2i(I, J, flavor=3, level=1, cut.min=cut.min, m=M, n=M, mat=FALSE)-> ind3
IND<- c( ind0, ind1, ind2, ind3)
NH<- length( IND)   # H0 will be NH by NH in size. 


# the Matern covariance function range=3, smoothness=1.
cov.obj<- matern.image.cov( grid=grid, setup=TRUE, theta=3, smoothness=1)

# find elements of D0= H0**2  using the formula
#  D0 = W Sigma W.t  where W is the  2-d W transform ( the inverse of W.i above).
#  and Sigma is the Matern covariance
#
# the following looping algorithms are a way to avoiding explicitly creating
# gigantic covariance matrices and wavelet basis matrices associated with large 
# grids. Of course in this example the grid is small enough (16X16) that the 
# matrices could be formed explicitly 
#

D0<- matrix( 0, NH,NH)
# fill in D0
 for ( k in 1:NH){
 tmp<- matrix( 0, M,M)
 tmp[IND[k]] <- 1
 hold<- Wtransform.image( tmp,  transpose=TRUE, cut.min=cut.min)
 hold<- matern.image.cov(Y=hold, cov.obj=cov.obj)
 hold<- Wtransform.image(hold, cut.min=cut.min)
 D0[k,] <- hold[IND]
}

# sqrt of D0
temp<- eigen( D0, symmetric=TRUE)

# next line should be H0<-temp$vectors%*%diag(sqrt(temp$values))%*%t(temp$vectors)

H0<- temp$vectors %*% diag(sqrt(temp$values))%*% t(temp$vectors)


# find H1
H1<- matrix(0,M,M)
for ( k in 1:MM){
 tmp<- matrix( 0, M,M)
 tmp[k] <- 1
 hold<- Wtransform.image( tmp,  transpose=TRUE, cut.min=cut.min)
 hold<- matern.image.cov(Y=hold, cov.obj=cov.obj)
 hold<- Wtransform.image(hold, cut.min=cut.min)
 H1[k] <- sqrt(hold[k]) }
# remember that H1 has the H0 diagonal values set to 1. 
H1[IND] <- 1

#OK good to go.  Create the H.obj list
H.obj<- list( m=M, n=M, ind0=IND, cut.min=cut.min, H0=H0, H1=H1)
#
#
#


# mutliply the covariance matrix by a random vector
tmp<- matrix( rnorm( 256), 16,16)
Wimage.cov( tmp, H.obj=H.obj)-> look

# generate a random field
Wimage.sim(H.obj)-> look
image.plot( look)


#A check that this really works!
#find  the  135 =(9-1)*16 +7 == (7,9) row of covariance matrix
# we know what this should look like. 

Wimage.cov( ind= 135, H.obj=H.obj, find.row=TRUE)-> look
image.plot( x,y,look) # the covariance of the grid point (7,9) with 
                      # all other grid points -- a bump centered at (7,9)

#multiply a vector by just a subset of Sigma

ind<- sample( 1:256,50) # random set of 50 indices
Y<- rnorm( 50) # random vector of length 50 
Wimage.cov(Y, ind= ind, H.obj=H.obj)[ind]-> look

# look is now of length 50 -- as expected. 
# In R notation, this finding Sigma[ind, ind] %*% Y

# OK suppose you really need Sigma[ind,ind] 
# e.g. in order for   solve( Sigma[ind,ind], u)  
# here is a loop to do it.

Sigma<- matrix( NA, 50,50)
for ( k in 1:50){
# for each member of subset find row and subset to just the columns identified by ind 
Sigma[k,] <- c( Wimage.cov( ind= ind[k], H.obj=H.obj, find.row=TRUE)[ind])
}