# exponential covariance matrix ( marginal variance =1) for the ozone
#locations
out<- Exp.cov( ChicagoO3$x, aRange=100)
# out is a 20X20 matrix
out2<- Exp.cov( ChicagoO3$x[6:20,],ChicagoO3$x[1:2,], aRange=100)
# out2 is 15X2 matrix
# Kriging fit where the nugget variance is found by GCV
# Matern covariance shape with range of 100.
#
fit<- Krig( ChicagoO3$x, ChicagoO3$y, Covariance="Matern", aRange=100,smoothness=2)
data( ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
# Omit the NAs
good<- !is.na( y)
x<- x[good,]
y<- y[good]
# example of calling the taper version directly
# Note that default covariance is exponential and default taper is
# Wendland (k=2).
stationary.taper.cov( x[1:3,],x[1:10,] , aRange=1.5, Taper.args= list(k=2,aRange=2.0,
dimension=2) )-> temp
# temp is now a tapered 3X10 cross covariance matrix in sparse format.
is.spam( temp) # evaluates to TRUE
# should be identical to
# the direct matrix product
temp2<- Exp.cov( x[1:3,],x[1:10,], aRange=1.5) * Wendland(rdist(x[1:3,],x[1:10,]),
aRange= 2.0, k=2, dimension=2)
test.for.zero( as.matrix(temp), temp2)
# Testing that the "V" option works as advertized ...
x1<- x[1:20,]
x2<- x[1:10,]
V<- matrix( c(2,1,0,4), 2,2)
Vi<- solve( V)
u1<- t(Vi%*% t(x1))
u2<- t(Vi%*% t(x2))
look<- exp(-1*rdist(u1,u2))
look2<- stationary.cov( x1,x2, V= V)
test.for.zero( look, look2)
# Here is an example of how the cross covariance multiply works
# and lots of options on the arguments
Ctest<- rnorm(10)
temp<- stationary.cov( x,x[1:10,], C= Ctest,
Covariance= "Wendland",
k=2, dimension=2, aRange=1.5 )
# do multiply explicitly
temp2<- stationary.cov( x,x[1:10,],
Covariance= "Wendland",
k=2, dimension=2, aRange=1.5 )%*% Ctest
test.for.zero( temp, temp2)
# use the tapered stationary version
# cov.args is part of the argument list passed to stationary.taper.cov
# within Krig.
# This example needs the spam package.
#
if (FALSE) {
Krig(x,y, cov.function = "stationary.taper.cov", aRange=1.5,
cov.args= list(Taper.args= list(k=2, dimension=2,aRange=2.0) )
) -> out2
# NOTE: Wendland is the default taper here.
}
# BTW this is very similar to
if (FALSE) {
Krig(x,y, aRange= 1.5)-> out
}
##################################################
#### nonstationary covariance Paciorek.cov
##################################################
if (FALSE) {
M<- 20
gridList<- list(x=seq( 0,1,length.out=M),
y=seq( 0,1,length.out=M))
sGrid<- make.surface.grid(gridList )
# An aRange surface
aRangeObj<- list(coef= c( 1,4,0))
class(aRangeObj)<- "myclass"
predict.myclass<- function( aRangeObj, x){
aRange<- exp(cbind( 1,x) %*% aRangeObj$coef)
return( aRange)
}
covMatrix<- Paciorek.cov( sGrid, sGrid, aRangeObj=aRangeObj)
# examine correlation surface between selected locations and the full grid.
set.panel( 2,2)
{
imagePlot( as.surface(sGrid, covMatrix[,10]))
imagePlot( as.surface(sGrid, covMatrix[,205]))
imagePlot( as.surface(sGrid, covMatrix[,305]))
imagePlot( as.surface(sGrid, covMatrix[,390]))
}
# simulation of the field
set.seed(222)
n<- nrow( sGrid)
f<- t(chol(covMatrix)) %*% rnorm(M^2)
set.panel()
imagePlot( as.surface(sGrid,f))
y<- f + .05*rnorm(n)
fitP<- spatialProcess( sGrid, y, cov.function="Paciorek.cov",
cov.args= list(aRangeObj = aRangeObj ) )
# check estimated tau and sigma
fitP$summary
# fitted surface
surface( fitP)
}
Run the code above in your browser using DataLab