R/bbase.bs.R

In DDPstar: Density Regression via Dirichlet Process Mixtures of Normal Structured Additive Regression Models

bbase.bs <-
function(x, ndx, bdeg = 3, pord = 2, eps = 1e-05, intercept = TRUE) {
	# B-spline basis
	xl <- min(x)
	xr <- max(x)
	dx <- (xr - xl)/ndx

	knots <- seq(xl - bdeg*dx, xr + bdeg*dx, by=dx)
	B <- splines::spline.des(knots, x, ord = bdeg + 1, outer.ok = TRUE)$design
	m <- ncol(B)

	# Penalty matrix
	Dp <- diff(diag(ncol(B)), differences = pord)

	# SVD of the penalty matrix	
	P.svd <- svd(crossprod(Dp))   # SVD penalty matrix  
	
	# Design matrix penalised part
	U.Z <- (P.svd$u)[,1:(m-pord)] # eigenvectors
	d <- (P.svd$d)[1:(m-pord)]    # eigenvalues
	Z <- B%*%U.Z

	# Design matrix unpenalised part
	X <- B%*%((P.svd$u)[,-(1:(m-pord))])
	D.temp <- sweep(X, 2, colMeans(X))				
	Xf <- svd(crossprod(D.temp))$u[,ncol(D.temp):1]
	X <- X%*%Xf
	U.X <- ((P.svd$u)[,-(1:(m-pord)), drop = FALSE])%*%Xf

	if(!intercept) {		
		# Remove intercept
		X.reduced <- X[,-1,drop = FALSE]
		
		# Design matrix
		B.new <- cbind(X.reduced, Z)
		# Center the matrix
		# cm <- colMeans(B.new)
		cm <- rep(0, ncol(B.new))
		B.new <- sweep(B.new, 2, cm)

		K <- list()
		K[[1]] <- diag(x = c(rep(0L, pord-1)), nrow = pord -1, ncol = pord -1)
		K[[2]] <- diag(x = d)	
	} else {	
		#U.X <- NULL
		#U.Z <- NULL
		#B.new <- B
		#K <- crossprod(Dp)
		#cm <- rep(0, ncol(B))
		
		# Design matrix
		B.new <- cbind(X,Z)
		# Center the matrix
		# cm <- colMeans(B.new)
		cm <- rep(0, ncol(B.new))
		B.new <- sweep(B.new, 2, cm)

		K <- list()
		K[[1]] <- diag(x = c(rep(0L, pord)), nrow = pord, ncol = pord)
		K[[2]] <- diag(x = d)

	}
	res <- list(B = B.new, B.orig = B, K = K, X = X, Z = Z, d = d)
	attr(res,"knots") <- knots
	attr(res,"bdeg") <- bdeg
	attr(res,"eps") <- eps
	attr(res,"intercept") <- intercept
	attr(res,"U.X") <- U.X
	attr(res,"U.Z") <- U.Z
	attr(res,"cm") <- cm
	class(res) <- c("bbase.bs", "matrix")
  	res
}