#we will use Franke's test function as an example f <- function(x, y) { .75*exp(-((9*x-2)^2 + (9*y-2)^2)/4) + .75*exp(-((9*x+1)^2/49 + (9*y+1)^2/10)) + .50*exp(-((9*x-7)^2 + (9*y-3)^2)/4) - .20*exp(-((9*x-4)^2 + (9*y-7)^2)) } #define a (fine) x-y grid and calculate the function values on the grid xf <- seq(1/26, 25/26, length=25); yf <- xf zf <- outer(xf, yf, f) #output all the plots that will be generated to the result.ps file (2 plots per page) #note that if you would like to take a look before output to the result.ps file, just #run the code between "begin plotting..." and "end of ...." postscript("result.ps", height=8, width=5, horizo=F) ## begin plotting the graphics par(mfrow=c(2,1)) #plot the Wendelberger's test function persp(xf, yf, zf, theta=130, phi=20, expand=0.45, xlab="x1", ylab="x2", zlab="f", ticktype="detailed", scale=F, main="true function") #generate a data set with Wendelberger's test function set.seed(223) N <- 13; xr <- (2*(1:N) - 1)/(2*N); yr <- xr zr <- outer(xr, yr, f); zrmax <- max(abs(zr)) #this is the noisy data we will use in the experiment, you may want to change 0.07 everywhere #and see what will happen noise <- rnorm(N^2, 0, 0.07*zrmax) zr <- zr + noise #plot the noisy data persp(xr, yr, zr, theta=130, phi=20, expand=0.45, xlab="x1", ylab="x2", zlab="y", xlim=c(0,1), ylim=c(0,1), zlim=range(zf), ticktype="detailed", scale=F, main="noisy data") #transpose the data to a column format to plug into 'Tps' function xc <- rep(xr, N) yc <- rep(yr, rep(N,N)) zc <- as.vector(zr) #or one could read the data from a data file #dat <- read.table("example.dat", header=T) #names(dat) <- c("xc", "yc", "zc") #attach(dat) #install and then load the fields package (containing Tps) and fit a thin plate spline #tuned by GCV install.packages("fields") library(fields) tpsfit <- Tps(cbind(xc, yc), zc) #predict the thin plate spline on the fine grid and plot the fitting ngrid <- length(xf) grid <- cbind(rep(xf, ngrid), rep(xf, rep(ngrid, ngrid))) out.p1 <- predict(tpsfit, grid) persp(xf, xf, matrix(out.p1, ngrid, ngrid, byrow=F), theta=130, phi=20, expand=0.45, xlab="x1", ylab="x2", zlab="y", xlim=c(0,1), ylim=c(0,1), zlim=range(zf), ticktype="detailed", scale=F, main="gcv fitting") #optional - look at lambda too big and lambda too small #use a bigger lambda, we specified lambda in "predict" function tpsfit2 <- Tps(cbind(xc, yc), zc) out.p2 <- predict(tpsfit2, grid,lambda=20*tpsfit$lambda) persp(xf, xf, matrix(out.p2, ngrid, ngrid, byrow=F), theta=130, phi=20, expand=0.45, xlab="x1", ylab="x2", zlab="y", xlim=c(0,1), ylim=c(0,1), zlim=range(zf), ticktype="detailed", scale=F, "oversmoothed") #use a smaller lambda , we specified lambda in "predict" function tpsfit3 <- Tps(cbind(xc, yc), zc) out.p3 <- predict(tpsfit3, grid, lambda=0.002*tpsfit$lambda) persp(xf, xf, matrix(out.p3, ngrid, ngrid, byrow=F), theta=130, phi=20, expand=0.45, xlab="x1", ylab="x2", zlab="y", xlim=c(0,1), ylim=c(0,1), zlim=range(zf), ticktype="detailed", scale=F, "undersmoothed") ## end of plotting the graphics #output to PS file graphics.off() #if you use motif(), you can shut down motif window by using dev.off() #output the smoothing parameter by GCV, the estimated variance, true variance #and the standard-deviation ratio (defined by estimated sd/true sd) #the ratio should be close to 1 name = c('lambda','estimated sigma2','true sigma2','sd-ratio') eq = c('=','=','=','=') val = c(tpsfit$lambda,tpsfit$shat.GCV^2,(0.07*zrmax)^2,tpsfit$shat.GCV/(0.07*zrmax)) out = data.frame(name,eq,val) write(t(out),file = 'result.txt',ncol = 3) #Plot the GCV curve postscript("gcv.ps", height=6, width=7, horizo=F) grid=tpsfit$gcv.grid plot(grid$GCV~log(grid$lambda),xlab='log(lambda)',ylab='GCV score', main='GCV curve',type='b') graphics.off() #if you hit tpsfit before quitting R with q(), you will #see other info sent by fields including degrees of freedom for signal and noise