Engauge Digitizer 2
Loading...
Searching...
No Matches
FittingStatistics.cpp
Go to the documentation of this file.
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"
9#include "FittingStatistics.h"
10#include "Logger.h"
11#include "Matrix.h"
12#include <QApplication>
13#include <qmath.h>
14
18
22
23void FittingStatistics::calculateCurveFit (int orderReduced,
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
37 while (!calculateCurveFitReducedFurther (orderReducedFurther,
39 significantDigits,
40 a)) {
41
42 // Retry if possible with lower order
44 if (orderReducedFurther < 0) {
45 break;
46 }
47 }
48 }
49
50 // Copy coefficients into member variable and into list for sending as a signal
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
69bool FittingStatistics::calculateCurveFitReducedFurther (int orderReducedFurther,
71 int significantDigits,
72 QVector<double> &a) const
73{
74 // Calculate X and y arrays in y = X a
77 loadXAndYArrays (orderReducedFurther,
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
86 LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFitReducedFurther determinant=" << denominator.determinant();
87
89 Matrix inv = denominator.inverse (significantDigits,
91
93
94 LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFitReducedFurther failed with order="
96
97 return false;
98 }
99
100 a = inv * X.transpose () * Y;
102
103 LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFitReducedFurther succeeded with order="
105 << " expectedIdentity="
106 << expectedIdentity.toString ().toLatin1().data ();
107
108 return true;
109}
110
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 (qFloor (order),
126 pointsConvenient.size() - 1);
127
128 calculateCurveFit (orderReduced,
131 significantDigits);
132 calculateStatistics (pointsConvenient,
134 mse,
135 rms,
136 rSquared);
137
138 qApp->restoreOverrideCursor();
139}
140
141void FittingStatistics::calculateStatistics (const FittingPointsConvenient &pointsConvenient,
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
175 }
176
177 mse = mseSum / pointsConvenient.count ();
178 rms = qSqrt (mse);
181 0);
182 }
183}
184
185void FittingStatistics::loadXAndYArrays (int orderReduced,
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
210double 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 ()) {
219 }
220 sum += coef * qPow (x, double (order));
221 }
222
223 return sum;
224}
const int INNER_RADIUS_MIN
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT if defined(QT_NO_DEBUG) && !defined(QT_FORCE_ASSERTS) define ENGAUGE...
QVector< double > FittingCurveCoefficients
Coefficients x0, x1, ... in y = a0 + a1 * x + a2 * x^2 + ...
QList< QPointF > FittingPointsConvenient
Array of (x,y) points in graph coordinates.
const int MAX_POLYNOMIAL_ORDER
log4cpp::Category * mainCat
Definition Logger.cpp:14
MatrixConsistent
Indicates if matrix is consistent (i.e. has at least one solution)
Definition Matrix.h:14
@ MATRIX_CONSISTENT
Definition Matrix.h:15
@ MATRIX_INCONSISTENT
Definition Matrix.h:16
FittingStatistics()
Single constructor.
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.
Matrix class that supports arbitrary NxN size.
Definition Matrix.h:21
double determinant() const
Return the determinant of this matrix.
Definition Matrix.cpp:66
Matrix inverse(int significantDigits, MatrixConsistent &matrixConsistent) const
Return the inverse of this matrix.
Definition Matrix.cpp:123
Matrix transpose() const
Return the transpose of the current matrix.
Definition Matrix.cpp:469
QString toString() const
Dump matrix to a string.
Definition Matrix.cpp:445
#define LOG4CPP_DEBUG_S(logger)
Definition convenience.h:20