LeastSquaresDegree2 <- function(t, b) { # Usage: t is independent variable vector, b is function data # i.e., b = f(t) # This code only works for finding second degree polynomial best-fit models t <- t(t) dataSize <- length(b) A <- mat.or.vec(dataSize, 3) # Built-in R function to create zero matrix or zero vector of arbitrary size # Given basis phi(z) = 1 + z + z^2 # Define matrix A A[1:dataSize,1] = 1 A[1:dataSize,2] = t A[1:dataSize,3] = t^2 # Compute QR decomposition of A. Pull Q and R out of QRdecomp QRdecomp <- qr(A) Q <- qr.Q(QRdecomp, complete=TRUE) R <- qr.R(QRdecomp, complete=TRUE) # Perform Q^* b^T (Conjugate transpose of Q) c <- Conj(t(Q))%*%t(b) # Find x. R isn't square - so we have to use qr.solve x <- qr.solve(R, c) cat("coefficients are:", x[1],",", x[2],",", x[3],".") # return(c(x=x)) # Now plot it. Least squares "l" plot first, then the points in red. plot(t, x[1] + x[2]*t + x[3]*t^2, type='l', xlab="independent variable t", ylab="function values f(t)", main="Data Plotted with 2nd Degree Least Squares Polynomial", col="blue") points(t, b, col="red") # For debugging # print(Q[190:201,190:201]) # print(dim(Q)) # print(dim(testQ)) # print(dim(b)) # testQ <- Conj(t(Q)) # print(testQ[190:201,190:201]) # print(c) # dim(c) # print(R) # Code by Eric Johnson -- Focus | Numeric }