[sysid] Fix adjusted R² calculation (#6101)

It hardcoded p to 2.

sysid/src/main/native/cpp/analysis/OLS.cpp

@@ -5,45 +5,84 @@
 #include "sysid/analysis/OLS.h"
 
 #include <cassert>
-#include <tuple>
-#include <vector>
+#include <cmath>
 
 #include <Eigen/Cholesky>
-#include <Eigen/Core>
 
-using namespace sysid;
+namespace sysid {
 
-std::tuple<std::vector<double>, double, double> sysid::OLS(
-    const Eigen::MatrixXd& X, const Eigen::VectorXd& y) {
+OLSResult OLS(const Eigen::MatrixXd& X, const Eigen::VectorXd& y) {
   assert(X.rows() == y.rows());
 
-  // The linear model can be written as follows:
-  // y = Xβ + u, where y is the dependent observed variable, X is the matrix
-  // of independent variables, β is a vector of coefficients, and u is a
-  // vector of residuals.
+  // The linear regression model can be written as follows:
+  //
+  //   y = Xβ + ε
+  //
+  // where y is the dependent observed variable, X is the matrix of independent
+  // variables, β is a vector of coefficients, and ε is a vector of residuals.
+  //
+  // We want to find the value of β that minimizes εᵀε.
+  //
+  //   ε = y − Xβ
+  //   εᵀε = (y − Xβ)ᵀ(y − Xβ)
+  //
+  //   β̂ = argmin (y − Xβ)ᵀ(y − Xβ)
+  //         β
+  //
+  // Take the partial derivative of the cost function with respect to β and set
+  // it equal to zero, then solve for β̂ .
+  //
+  //   0 = −2Xᵀ(y − Xβ̂)
+  //   0 = Xᵀ(y − Xβ̂)
+  //   0 = Xᵀy − XᵀXβ̂
+  //   XᵀXβ̂ = Xᵀy
+  //   β̂  = (XᵀX)⁻¹Xᵀy
 
-  // We want to minimize u² = uᵀu = (y - Xβ)ᵀ(y - Xβ).
   // β = (XᵀX)⁻¹Xᵀy
+  //
+  // XᵀX is guaranteed to be symmetric positive definite, so an LLT
+  // decomposition can be used.
+  Eigen::MatrixXd β = (X.transpose() * X).llt().solve(X.transpose() * y);
 
-  // Calculate β that minimizes uᵀu.
-  Eigen::MatrixXd beta = (X.transpose() * X).llt().solve(X.transpose() * y);
+  // Error sum of squares
+  double SSE = (y - X * β).squaredNorm();
 
-  // We will now calculate R² or the coefficient of determination, which
-  // tells us how much of the total variation (variation in y) can be
-  // explained by the regression model.
+  // Sample size
+  int n = X.rows();
 
-  // We will first calculate the sum of the squares of the error, or the
-  // variation in error (SSE).
-  double SSE = (y - X * beta).squaredNorm();
+  // Number of explanatory variables
+  int p = β.rows();
 
-  int n = X.cols();
+  // Total sum of squares (total variation in y)
+  //
+  // From slide 24 of
+  // http://www.stat.columbia.edu/~fwood/Teaching/w4315/Fall2009/lecture_11:
+  //
+  //   SSTO = yᵀy - 1/n yᵀJy
+  //
+  // Let J = I.
+  //
+  //   SSTO = yᵀy - 1/n yᵀy
+  //   SSTO = (n − 1)/n yᵀy
+  double SSTO = (n - 1.0) / n * (y.transpose() * y).value();
 
-  // Now we will calculate the total variation in y, known as SSTO.
-  double SSTO = ((y.transpose() * y) - (1.0 / n) * (y.transpose() * y)).value();
+  // R² or the coefficient of determination, which represents how much of the
+  // total variation (variation in y) can be explained by the regression model
+  double rSquared = 1.0 - SSE / SSTO;
 
-  double rSquared = (SSTO - SSE) / SSTO;
-  double adjRSquared = 1.0 - (1.0 - rSquared) * ((n - 1.0) / (n - 3.0));
+  // Adjusted R²
+  //
+  //                       n − 1
+  //   R̅² = 1 − (1 − R²) ---------
+  //                     n − p − 1
+  //
+  // See https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2
+  double adjRSquared = 1.0 - (1.0 - rSquared) * ((n - 1.0) / (n - p - 1.0));
+
+  // Root-mean-square error
   double RMSE = std::sqrt(SSE / n);
 
-  return {{beta.data(), beta.data() + beta.rows()}, adjRSquared, RMSE};
+  return {{β.data(), β.data() + β.size()}, adjRSquared, RMSE};
 }
+
+}  // namespace sysid