refund (version 0.1-37)

fbps: Sandwich smoother for matrix data

Description

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

Usage

fbps(
  data,
  subj = NULL,
  covariates = NULL,
  knots = 35,
  knots.option = "equally-spaced",
  periodicity = c(FALSE, FALSE),
  p = 3,
  m = 2,
  lambda = NULL,
  selection = "GCV",
  search.grid = T,
  search.length = 100,
  method = "L-BFGS-B",
  lower = -20,
  upper = 20,
  control = NULL
)

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

Arguments

data

n1 by n2 data matrix without missing data

subj

vector of subject id (corresponding to the columns of data); defaults to NULL

covariates

list of two vectors of covariates of lengths 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

knots.option

knot selection method; defaults to "equally-spaced"

periodicity

vector of two logical, indicating periodicity in the direction of row and column; defaults to c(FALSE, FALSE)

p

degrees of B-splines; defaults to 3

m

order of differencing penalty; defaults to 2

lambda

user-specified smoothing parameters; defaults to NULL

selection

selection of smoothing parameter; defaults to "GCV"

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

Author

Luo Xiao

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,"\n")
cat("The smoothing parameters are:",est$lambda,"\n")
########################################################################
########## 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")