00001
00002
00003
00004
00005
00006
00007 #include "EngaugeAssert.h"
00008 #include "FittingCurveCoefficients.h"
00009 #include "FittingStatistics.h"
00010 #include "Logger.h"
00011 #include "Matrix.h"
00012 #include <QApplication>
00013 #include <qmath.h>
00014
00015 FittingStatistics::FittingStatistics ()
00016 {
00017 }
00018
00019 FittingStatistics::~FittingStatistics()
00020 {
00021 }
00022
00023 void FittingStatistics::calculateCurveFit (int orderReduced,
00024 const FittingPointsConvenient &pointsConvenient,
00025 FittingCurveCoefficients &coefficients)
00026 {
00027 QVector<double> a;
00028
00029 if (0 <= orderReduced) {
00030
00031
00032 Matrix X (pointsConvenient.size (), orderReduced + 1);
00033 QVector<double> Y (pointsConvenient.size ());
00034 loadXAndYArrays (orderReduced,
00035 pointsConvenient,
00036 X,
00037 Y);
00038
00039
00040 Matrix denominator = X.transpose () * X;
00041 LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFit determinant=" << denominator.determinant();
00042 a = denominator.inverse () * X.transpose () * Y;
00043
00044 Matrix expectedIdentity = denominator * denominator.inverse ();
00045 LOG4CPP_DEBUG_S ((*mainCat)) << "FittingStatistics::calculateCurveFit expectedIdentity="
00046 << expectedIdentity.toString ().toLatin1().data ();
00047 }
00048
00049
00050 FittingCurveCoefficients fittingCurveCoef;
00051 for (int order = 0; order <= MAX_POLYNOMIAL_ORDER; order++) {
00052
00053 if (order <= orderReduced) {
00054
00055
00056 coefficients [order] = a [order];
00057 fittingCurveCoef.append (a [order]);
00058
00059 } else {
00060
00061
00062 coefficients [order] = 0;
00063
00064 }
00065 }
00066 }
00067
00068 void FittingStatistics::calculateCurveFitAndStatistics (unsigned int order,
00069 const FittingPointsConvenient &pointsConvenient,
00070 FittingCurveCoefficients &coefficients,
00071 double &mse,
00072 double &rms,
00073 double &rSquared)
00074 {
00075
00076 qApp->setOverrideCursor (Qt::WaitCursor);
00077
00078
00079
00080
00081 int orderReduced = qMin ((int) order,
00082 pointsConvenient.size() - 1);
00083
00084 calculateCurveFit (orderReduced,
00085 pointsConvenient,
00086 coefficients);
00087 calculateStatistics (pointsConvenient,
00088 coefficients,
00089 mse,
00090 rms,
00091 rSquared);
00092
00093 qApp->restoreOverrideCursor();
00094 }
00095
00096 void FittingStatistics::calculateStatistics (const FittingPointsConvenient &pointsConvenient,
00097 const FittingCurveCoefficients &coefficients,
00098 double &mse,
00099 double &rms,
00100 double &rSquared)
00101 {
00102
00103 double ySum = 0;
00104 FittingPointsConvenient::const_iterator itrC;
00105 for (itrC = pointsConvenient.begin (); itrC != pointsConvenient.end (); itrC++) {
00106
00107 const QPointF &pointC = *itrC;
00108 ySum += pointC.y();
00109 }
00110 double yAverage = ySum / pointsConvenient.length();
00111
00112
00113 mse = 0;
00114 rms = 0;
00115 rSquared = 0;
00116
00117 if (pointsConvenient.count() > 0) {
00118
00119 double mseSum = 0, rSquaredNumerator = 0, rSquaredDenominator = 0;
00120 for (itrC = pointsConvenient.begin(); itrC != pointsConvenient.end (); itrC++) {
00121
00122 const QPointF &pointC = *itrC;
00123 double yActual = pointC.y();
00124 double yCurveFit = yFromXAndCoefficients (coefficients,
00125 pointC.x());
00126
00127 mseSum += (yCurveFit - yActual ) * (yCurveFit - yActual );
00128 rSquaredNumerator += (yCurveFit - yAverage) * (yCurveFit - yAverage);
00129 rSquaredDenominator += (yActual - yAverage) * (yActual - yAverage);
00130 }
00131
00132 mse = mseSum / pointsConvenient.count ();
00133 rms = qSqrt (mse);
00134 rSquared = (rSquaredDenominator > 0 ?
00135 rSquaredNumerator / rSquaredDenominator :
00136 0);
00137 }
00138 }
00139
00140 void FittingStatistics::loadXAndYArrays (int orderReduced,
00141 const FittingPointsConvenient &pointsConvenient,
00142 Matrix &X,
00143 QVector<double> &Y) const
00144 {
00145 ENGAUGE_ASSERT (Y.size () == X.rows ());
00146
00147
00148 int row;
00149 FittingPointsConvenient::const_iterator itr;
00150 for (row = 0, itr = pointsConvenient.begin(); itr != pointsConvenient.end(); itr++, row++) {
00151
00152 const QPointF &p = *itr;
00153 double x = p.x ();
00154 double y = p.y ();
00155
00156 for (int order = 0; order <= orderReduced; order++) {
00157
00158 X.set (row, order, qPow (x, order));
00159 }
00160
00161 Y [row] = y;
00162 }
00163 }
00164
00165 double FittingStatistics::yFromXAndCoefficients (const FittingCurveCoefficients &coefficients,
00166 double x) const
00167 {
00168 double sum = 0;
00169
00170 for (int order = 0; order <= MAX_POLYNOMIAL_ORDER; order++) {
00171 sum += coefficients [order] * qPow (x, (double) order);
00172 }
00173
00174 return sum;
00175 }