Engauge Digitizer  2
 All Classes Functions Variables Typedefs Enumerations Friends Pages
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  int significantDigits)
27 {
28  QVector<double> a; // Unused if there are no points
29 
30  if (0 <= orderReduced) {
31 
32  // We further reduce the order in case the inputs are badly defined. Known bad data example
33  // is trying to linearly fit two points with same x but different y values (=infinite slope).
34  // Retries are done with successively reduced order. With good inputs no reduction is required
35  // so orderReducedFurther succeeds with the same value as orderReduced
36  int orderReducedFurther = orderReduced;
37  while (!calculateCurveFitReducedFurther (orderReducedFurther,
38  pointsConvenient,
39  significantDigits,
40  a)) {
41 
42  // Retry if possible with lower order
43  --orderReducedFurther;
44  if (orderReducedFurther < 0) {
45  break;
46  }
47  }
48  }
49 
50  // Copy coefficients into member variable and into list for sending as a signal
51  FittingCurveCoefficients fittingCurveCoef;
52  for (int order = 0; order <= MAX_POLYNOMIAL_ORDER; order++) {
53 
54  if (order < a.size ()) {
55 
56  // Copy from polynomial regression vector
57  coefficients [order] = a [order];
58  fittingCurveCoef.append (a [order]);
59 
60  } else if (order < coefficients.size ()) {
61 
62  // Set to zero in case value gets used somewhere
63  coefficients [order] = 0;
64 
65  }
66  }
67 }
68 
69 bool FittingStatistics::calculateCurveFitReducedFurther (int orderReducedFurther,
70  const FittingPointsConvenient &pointsConvenient,
71  int significantDigits,
72  QVector<double> &a) const
73 {
74  // Calculate X and y arrays in y = X a
75  Matrix X (pointsConvenient.size (), orderReducedFurther + 1);
76  QVector<double> Y (pointsConvenient.size ());
77  loadXAndYArrays (orderReducedFurther,
78  pointsConvenient,
79  X,
80  Y);
81 
82  // Solve for the coefficients a in y = X a + epsilon using a = (Xtranpose X)^(-1) Xtranspose y. Solution
83  // uses more complicated transposes rather than just using a = X^(-1) y since the transpose approach
84  // can handle when the a matrix is nonsquare
85  Matrix denominator = X.transpose () * X;
86  LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFitReducedFurther determinant=" << denominator.determinant();
87 
88  MatrixConsistent matrixConsistent = MATRIX_CONSISTENT;
89  Matrix inv = denominator.inverse (significantDigits,
90  matrixConsistent);
91 
92  if (matrixConsistent == MATRIX_INCONSISTENT) {
93 
94  LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFitReducedFurther failed with order="
95  << orderReducedFurther;
96 
97  return false;
98  }
99 
100  a = inv * X.transpose () * Y;
101  Matrix expectedIdentity = denominator * inv;
102 
103  LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFitReducedFurther succeeded with order="
104  << orderReducedFurther
105  << " expectedIdentity="
106  << expectedIdentity.toString ().toLatin1().data ();
107 
108  return true;
109 }
110 
112  const FittingPointsConvenient &pointsConvenient,
113  FittingCurveCoefficients &coefficients,
114  double &mse,
115  double &rms,
116  double &rSquared,
117  int significantDigits)
118 {
119  // Let user know something is happening if a high order was picked since that can take a long time
120  qApp->setOverrideCursor (Qt::WaitCursor);
121 
122  // To prevent having an underdetermined system with an infinite number of solutions (which will result
123  // in divide by zero when computing an inverse) we reduce the order here if necessary.
124  // In other words, we limit the order to -1 for no points, 0 for one point, 1 for two points, and so on
125  int orderReduced = qMin ((int) order,
126  pointsConvenient.size() - 1);
127 
128  calculateCurveFit (orderReduced,
129  pointsConvenient,
130  coefficients,
131  significantDigits);
132  calculateStatistics (pointsConvenient,
133  coefficients,
134  mse,
135  rms,
136  rSquared);
137 
138  qApp->restoreOverrideCursor();
139 }
140 
141 void FittingStatistics::calculateStatistics (const FittingPointsConvenient &pointsConvenient,
142  const FittingCurveCoefficients &coefficients,
143  double &mse,
144  double &rms,
145  double &rSquared)
146 {
147  // First pass to get average y
148  double ySum = 0;
149  FittingPointsConvenient::const_iterator itrC;
150  for (itrC = pointsConvenient.begin (); itrC != pointsConvenient.end (); itrC++) {
151 
152  const QPointF &pointC = *itrC;
153  ySum += pointC.y();
154  }
155  double yAverage = ySum / pointsConvenient.length();
156 
157  // Second pass to compute squared terms
158  mse = 0;
159  rms = 0;
160  rSquared = 0;
161 
162  if (pointsConvenient.count() > 0) {
163 
164  double mseSum = 0, rSquaredNumerator = 0, rSquaredDenominator = 0;
165  for (itrC = pointsConvenient.begin(); itrC != pointsConvenient.end (); itrC++) {
166 
167  const QPointF &pointC = *itrC;
168  double yActual = pointC.y();
169  double yCurveFit = yFromXAndCoefficients (coefficients,
170  pointC.x());
171 
172  mseSum += (yCurveFit - yActual ) * (yCurveFit - yActual );
173  rSquaredNumerator += (yCurveFit - yAverage) * (yCurveFit - yAverage);
174  rSquaredDenominator += (yActual - yAverage) * (yActual - yAverage);
175  }
176 
177  mse = mseSum / pointsConvenient.count ();
178  rms = qSqrt (mse);
179  rSquared = (rSquaredDenominator > 0 ?
180  rSquaredNumerator / rSquaredDenominator :
181  0);
182  }
183 }
184 
185 void FittingStatistics::loadXAndYArrays (int orderReduced,
186  const FittingPointsConvenient &pointsConvenient,
187  Matrix &X,
188  QVector<double> &Y) const
189 {
190  ENGAUGE_ASSERT (Y.size () == X.rows ());
191 
192  // Construct the matrices
193  int row;
194  FittingPointsConvenient::const_iterator itr;
195  for (row = 0, itr = pointsConvenient.begin(); itr != pointsConvenient.end(); itr++, row++) {
196 
197  const QPointF &p = *itr;
198  double x = p.x ();
199  double y = p.y ();
200 
201  for (int order = 0; order <= orderReduced; order++) {
202 
203  X.set (row, order, qPow (x, order));
204  }
205 
206  Y [row] = y;
207  }
208 }
209 
210 double FittingStatistics::yFromXAndCoefficients (const FittingCurveCoefficients &coefficients,
211  double x) const
212 {
213  double sum = 0;
214 
215  for (int order = 0; order <= MAX_POLYNOMIAL_ORDER; order++) {
216  double coef = 0;
217  if (order < coefficients.size ()) {
218  coef = coefficients [order];
219  }
220  sum += coef * qPow (x, (double) order);
221  }
222 
223  return sum;
224 }
int rows() const
Height of matrix.
Definition: Matrix.cpp:423
Matrix inverse(int significantDigits, MatrixConsistent &matrixConsistent) const
Return the inverse of this matrix.
Definition: Matrix.cpp:123
double determinant() const
Return the determinant of this matrix.
Definition: Matrix.cpp:66
void calculateCurveFitAndStatistics(unsigned int order, const FittingPointsConvenient &pointsConvenient, FittingCurveCoefficients &coefficients, double &mse, double &rms, double &rSquared, int significantDigits)
Compute the curve fit and the statistics for that curve fit.
FittingStatistics()
Single constructor.
Matrix transpose() const
Return the transpose of the current matrix.
Definition: Matrix.cpp:469
Matrix class that supports arbitrary NxN size.
Definition: Matrix.h:20
QString toString() const
Dump matrix to a string.
Definition: Matrix.cpp:445
void set(int row, int col, double value)
Set (row, col) element.
Definition: Matrix.cpp:428