refund (version 0.1-1)

fbps: Sandwich smoother for matrix data

Description

A fast bivariate P-spline method for smoothing matrix data.

Usage

fbps(data, covariates = NULL, knots = 35, p = 3, m = 2, lambda = NULL, 
     alpha = 1, search.grid = T, search.length = 100, method = "L-BFGS-B",
     lower = -20, upper = 20, control = NULL)

Arguments

data

n1 by n2 data matrix without missing data

covariates

list of two vectors of covariates of lengthes n1 and n2; if NULL, then generates equidistant covariates

knots

list of two vectors of knots or number of equidistant knots for all dimensions; defaults to 35

p

degrees of B-splines; defaults to 3

m

order of differencing penalty; defaults to 2

lambda

user-specified smoothing parameters; defaults to NULL

alpha

the coefficient of the penalty for degrees of freedom in the GCV criterion; defaults to 1

search.grid

logical; defaults to TRUE, if FALSE, uses optim

search.length

number of equidistant (log scale) smoothing parameter; defaults to 100

method

see optim; defaults to L-BFGS-B

lower, upper

bounds for log smoothing parameter, passed to optim; defaults are -20 and 20.

control

see optim

Value

• A list with components
• lambda
  vector of length 2 of selected smoothing parameters
• Yhat
  fitted data
• trace
  trace of the overall smoothing matrix
• gcv
  value of generalized cross validation
• Theta
  matrix of estimated coefficients

Details

The smoothing parameter can be user-specified; otherwise, the function uses grid searching method or optim for selecting the smoothing parameter.

References

Xiao, L., Li, Y., and Ruppert, D. (2013). Fast bivariate P-splines: the sandwich smoother. Journal of the Royal Statistical Society: Series B, 75(3), 577--599.

Examples

##########################
#### True function   #####
##########################
n1 <- 60
n2 <- 80
x <- (1: n1)/n1-1/2/n1
z <- (1: n2)/n2-1/2/n2
MY <- array(0,c(length(x),length(z)))

sigx <- .3
sigz <- .4
for(i in 1: length(x))
for(j in 1: length(z))
{
#MY[i,j] <- .75/(pi*sigx*sigz) *exp(-(x[i]-.2)^2/sigx^2-(z[j]-.3)^2/sigz^2)
#MY[i,j] <- MY[i,j] + .45/(pi*sigx*sigz) *exp(-(x[i]-.7)^2/sigx^2-(z[j]-.8)^2/sigz^2)
MY[i,j] = sin(2*pi*(x[i]-.5)^3)*cos(4*pi*z[j])
}
##########################
#### Observed data   #####
##########################
sigma <- 1
Y <- MY + sigma*rnorm(n1*n2,0,1)
##########################
####   Estimation    #####
##########################

est <- fbps(Y,list(x=x,z=z))
mse <- mean((est$Yhat-MY)^2)
cat("mse of fbps is",mse,"")
cat("The smoothing parameters are:",est$lambda,"")
########################################################################
########## Compare the estimated surface with the true surface #########
########################################################################

par(mfrow=c(1,2))
persp(x,z,MY,zlab="f(x,z)",zlim=c(-1,2.5), phi=30,theta=45,expand=0.8,r=4,
      col="blue",main="True surface")
persp(x,z,est$Yhat,zlab="f(x,z)",zlim=c(-1,2.5),phi=30,theta=45,
      expand=0.8,r=4,col="red",main="Estimated surface")