R/trapzmat.R

In fda: Functional Data Analysis

trapzmat <- function(X,Y,delta=1,wt=rep(1,n)) {
#TRAPZMAT integrates the products of two matrices of values
#   using the trapezoidal rule, assuming equal spacing
#  X is the first  matrix of values
#  Y is the second matrix of values
#  DELTA is the spacing between argument values (one by default)
#  WT is a vector of weights (ones by default)
#
#  XtWY is a matrix of integral estimates, number of rows equal to
#  number of col of X, number of cols equal to number of cols of Y

	X <- as.matrix(X)
	Y <- as.matrix(Y)
	
	n <- dim(X)[1]

	if (dim(Y)[1] != n) {
    	stop("X and Y do not have same number of rows.")
	}

	if (length(wt) != n) {
    	stop("X and WT do not have same number of rows.")
	}

	if (delta <= 0) {
    	stop("DELTA is not a positive value.")
	}

	wt[c(1,n)] <- wt[c(1,n)]/2
	wt <- wt*delta

	X <- X*outer(wt,rep(1,dim(X)[2]))
	XtWY <- crossprod(X,Y)
	return(XtWY)
}