fbps | R Documentation |
A fast bivariate P-spline method for smoothing matrix data.
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
)
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
|
search.length |
number of equidistant (log scale) smoothing parameter; defaults to 100 |
method |
see |
lower , upper |
bounds for log smoothing parameter, passed to
|
control |
see |
The smoothing parameter can be user-specified; otherwise, the function uses
grid searching method or optim
for selecting the smoothing
parameter.
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 |
Luo Xiao lxiao@jhsph.edu
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.
##########################
#### 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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.