From 2a6083071f3c3a3757b8136f702d944c58a8a67a Mon Sep 17 00:00:00 2001
From: Mark Liversedge <liversedge@gmail.com>
Date: Sat, 25 Aug 2018 07:47:23 +0100
Subject: [PATCH] LR Model Fit calculate R2

.. lots of folks understand what R2 means for a linear
   regression so calculate it alongside RMSE when we
   are performing a linear regression.
---
 src/Charts/CPPlot.cpp   |  2 +-
 src/Metrics/PDModel.cpp | 23 +++++++++++++++++++----
 2 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/Charts/CPPlot.cpp b/src/Charts/CPPlot.cpp
index 7a9ad7e69..43556b288 100644
--- a/src/Charts/CPPlot.cpp
+++ b/src/Charts/CPPlot.cpp
@@ -703,7 +703,7 @@ CPPlot::plotModel()
     // Pmax is often higher than the test values (they're for
     // 3-20 mins typically so well short of pmax).
     if (!showDelta && rideSeries == RideFile::watts && pdModel && pdModel->PMax() > ymax) {
-        if (pdModel->PMax() > ymax) setAxisScale(yLeft, 0, pdModel->PMax());
+        if (pdModel->PMax() > ymax) setAxisScale(yLeft, 0, pdModel->PMax() * 1.1f);
     }
 }
 
diff --git a/src/Metrics/PDModel.cpp b/src/Metrics/PDModel.cpp
index a66ebf2f8..bac2dbe47 100644
--- a/src/Metrics/PDModel.cpp
+++ b/src/Metrics/PDModel.cpp
@@ -176,11 +176,26 @@ PDModel::deriveCPParameters(int model)
         // get vector of residuals
         QVector<double> residuals(t.size());
         double errtot=0;
+        double sse=0;
+        double meany=0;
         for(int i=0; i<t.size(); i++){
-            double err = p[i] - y(t[i]/60.0f);
-            errtot += err; // for mean
-            residuals[i]=err;
+            double py = y(t[i]/60.0f);  // predicted y
+            meany += p[i];              // calculatng mean divided at end
+            double err = p[i] - py;     // error for this point
+            errtot += err;              // for mean errors divided at end
+            sse += err*err;             // sum of squared errors for r2 calc
+            residuals[i]=err;           // vector of residuals
         }
+
+        // lets use to calc R2
+        meany /= double(t.size());
+        double sv=0;
+        for(int i=0; i<t.size(); i++) {
+            sv += (p[i]-meany) * (p[i]-meany);
+        }
+        double R2=1-(sse/sv);
+
+        // carry on with RMSE
         double mean = errtot/double(t.size());
         errtot = 0;
 
@@ -190,7 +205,7 @@ PDModel::deriveCPParameters(int model)
 
         // RMSE
         double RMSE=sqrt(mean);
-        fitsummary = QString("RMSE %1 watts [LR] %2 points").arg(RMSE, 0, 'f', 2).arg(t.size());
+        fitsummary = QString("RMSE %1 watts R2=%3 [LR] %2 points").arg(RMSE, 0, 'f', 2).arg(t.size()).arg(R2, 0, 'f', 3);
 
     } else if (fit == LeastSquares && this->nparms() > 0) {
     // only try lmfit if the model supports that