Engauge Digitizer  2
FittingStatistics.cpp
1 /******************************************************************************************************
2  * (C) 2016 markummitchell@github.com. This file is part of Engauge Digitizer, which is released *
3  * under GNU General Public License version 2 (GPLv2) or (at your option) any later version. See file *
4  * LICENSE or go to gnu.org/licenses for details. Distribution requires prior written permission. *
5  ******************************************************************************************************/
6 
7 #include "EngaugeAssert.h"
8 #include "FittingCurveCoefficients.h"
9 #include "FittingStatistics.h"
10 #include "Logger.h"
11 #include "Matrix.h"
12 #include <QApplication>
13 #include <qmath.h>
14 
16 {
17 }
18 
19 FittingStatistics::~FittingStatistics()
20 {
21 }
22 
23 void FittingStatistics::calculateCurveFit (int orderReduced,
24  const FittingPointsConvenient &pointsConvenient,
25  FittingCurveCoefficients &coefficients)
26 {
27  QVector<double> a; // Unused if there are no points
28 
29  if (0 <= orderReduced) {
30 
31  // Calculate X and y arrays in y = X a
32  Matrix X (pointsConvenient.size (), orderReduced + 1);
33  QVector<double> Y (pointsConvenient.size ());
34  loadXAndYArrays (orderReduced,
35  pointsConvenient,
36  X,
37  Y);
38 
39  // Solve for the coefficients a in y = X a + epsilon using a = (Xtranpose X)^(-1) Xtranspose y
40  Matrix denominator = X.transpose () * X;
41  LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFit determinant=" << denominator.determinant();
42  a = denominator.inverse () * X.transpose () * Y;
43 
44  Matrix expectedIdentity = denominator * denominator.inverse ();
45  LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFit expectedIdentity="
46  << expectedIdentity.toString ().toLatin1().data ();
47  }
48 
49  // Copy coefficients into member variable and into list for sending as a signal
50  FittingCurveCoefficients fittingCurveCoef;
51  for (int order = 0; order <= MAX_POLYNOMIAL_ORDER; order++) {
52 
53  if (order <= orderReduced) {
54 
55  // Copy from polynomial regression vector
56  coefficients [order] = a [order];
57  fittingCurveCoef.append (a [order]);
58 
59  } else {
60 
61  // Set to zero in case value gets used somewhere
62  coefficients [order] = 0;
63 
64  }
65  }
66 }
67 
69  const FittingPointsConvenient &pointsConvenient,
70  FittingCurveCoefficients &coefficients,
71  double &mse,
72  double &rms,
73  double &rSquared)
74 {
75  // Let user know something is happening if a high order was picked since that can take a long time
76  qApp->setOverrideCursor (Qt::WaitCursor);
77 
78  // To prevent having an underdetermined system with an infinite number of solutions (which will result
79  // in divide by zero when computing an inverse) we reduce the order here if necessary.
80  // In other words, we limit the order to -1 for no points, 0 for one point, 1 for two points, and so on
81  int orderReduced = qMin ((int) order,
82  pointsConvenient.size() - 1);
83 
84  calculateCurveFit (orderReduced,
85  pointsConvenient,
86  coefficients);
87  calculateStatistics (pointsConvenient,
88  coefficients,
89  mse,
90  rms,
91  rSquared);
92 
93  qApp->restoreOverrideCursor();
94 }
95 
96 void FittingStatistics::calculateStatistics (const FittingPointsConvenient &pointsConvenient,
97  const FittingCurveCoefficients &coefficients,
98  double &mse,
99  double &rms,
100  double &rSquared)
101 {
102  // First pass to get average y
103  double ySum = 0;
104  FittingPointsConvenient::const_iterator itrC;
105  for (itrC = pointsConvenient.begin (); itrC != pointsConvenient.end (); itrC++) {
106 
107  const QPointF &pointC = *itrC;
108  ySum += pointC.y();
109  }
110  double yAverage = ySum / pointsConvenient.length();
111 
112  // Second pass to compute squared terms
113  mse = 0;
114  rms = 0;
115  rSquared = 0;
116 
117  if (pointsConvenient.count() > 0) {
118 
119  double mseSum = 0, rSquaredNumerator = 0, rSquaredDenominator = 0;
120  for (itrC = pointsConvenient.begin(); itrC != pointsConvenient.end (); itrC++) {
121 
122  const QPointF &pointC = *itrC;
123  double yActual = pointC.y();
124  double yCurveFit = yFromXAndCoefficients (coefficients,
125  pointC.x());
126 
127  mseSum += (yCurveFit - yActual ) * (yCurveFit - yActual );
128  rSquaredNumerator += (yCurveFit - yAverage) * (yCurveFit - yAverage);
129  rSquaredDenominator += (yActual - yAverage) * (yActual - yAverage);
130  }
131 
132  mse = mseSum / pointsConvenient.count ();
133  rms = qSqrt (mse);
134  rSquared = (rSquaredDenominator > 0 ?
135  rSquaredNumerator / rSquaredDenominator :
136  0);
137  }
138 }
139 
140 void FittingStatistics::loadXAndYArrays (int orderReduced,
141  const FittingPointsConvenient &pointsConvenient,
142  Matrix &X,
143  QVector<double> &Y) const
144 {
145  ENGAUGE_ASSERT (Y.size () == X.rows ());
146 
147  // Construct the matrices
148  int row;
149  FittingPointsConvenient::const_iterator itr;
150  for (row = 0, itr = pointsConvenient.begin(); itr != pointsConvenient.end(); itr++, row++) {
151 
152  const QPointF &p = *itr;
153  double x = p.x ();
154  double y = p.y ();
155 
156  for (int order = 0; order <= orderReduced; order++) {
157 
158  X.set (row, order, qPow (x, order));
159  }
160 
161  Y [row] = y;
162  }
163 }
164 
165 double FittingStatistics::yFromXAndCoefficients (const FittingCurveCoefficients &coefficients,
166  double x) const
167 {
168  double sum = 0;
169 
170  for (int order = 0; order <= MAX_POLYNOMIAL_ORDER; order++) {
171  sum += coefficients [order] * qPow (x, (double) order);
172  }
173 
174  return sum;
175 }
Matrix inverse() const
Return the inverse of this matrix.
Definition: Matrix.cpp:121
Matrix transpose() const
Return the transpose of the current matrix.
Definition: Matrix.cpp:410
double determinant() const
Return the determinant of this matrix.
Definition: Matrix.cpp:65
FittingStatistics()
Single constructor.
Matrix class that supports arbitrary NxN size.
Definition: Matrix.h:14
int rows() const
Height of matrix.
Definition: Matrix.cpp:364
void set(int row, int col, double value)
Set (row, col) element.
Definition: Matrix.cpp:369
void calculateCurveFitAndStatistics(unsigned int order, const FittingPointsConvenient &pointsConvenient, FittingCurveCoefficients &coefficients, double &mse, double &rms, double &rSquared)
Compute the curve fit and the statistics for that curve fit.
QString toString() const
Dump matrix to a string.
Definition: Matrix.cpp:386