Changeset 12338

Timestamp:

Mar 8, 2007, 1:23:00 PM (19 years ago)

Author:

rhl

Message:

Allow for loss of degrees of freedom in fitting when estimating variance

File:

• branches/rel-1_0/psLib/src/math/psMinimizePolyFit.c (modified) (6 diffs)

branches/rel-1_0/psLib/src/math/psMinimizePolyFit.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-01-09 22:38:53 $
+ *  @version $Revision: 1.29.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-03-08 23:23:00 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -879 +879 @@
         double stdevValue = psStatsGetValue (stats, stdevOption);
 
+        // correct for the number of degrees of freedom absorbed by the fit
+        const int nterm = psPolynomial1DNterm(poly);
+        if (Nkeep <= nterm) {   // we can't estimate the standard deviation
+            psError(PS_ERR_UNKNOWN, true, "Number of data points %d is <= number of degrees of freedom %d",
+                    Nkeep, nterm);
+            psFree(resid);
+            psFree(fit);
+            
+            return false;
+        }
+        stdevValue *= (Nkeep - 1)/(Nkeep - nterm - 1);
+
         psTrace("psLib.math", 5, "Mean is %f\n", meanValue);
         psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue);
@@ -1045 +1057 @@
     }
     psFree(xySums);
-
+#if 0
+    /*
+     * psMatrixGJSolve doesn't detect singular matrices very well. I'll do an SVN decomp
+     * to test for this -- but this is quite expensive, and a better solver would be wise. RHL.
+     *
+     * Don't run this test for now as I made the Gauss-Jordan solver require that |pivot| > FLT_EPSILON
+     *
+     * We could use lambda/evectors to solve the system without proceeding to a Gauss-Jordan
+     * solver, but I'm too lazy for now
+     */
+    psVector *lambda = psVectorAlloc(A->numCols, A->type.type); // eigenvalues
+    psImage *evectors = psMatrixSVD(NULL, lambda, A); // eigenvectors
+
+    for (int i = 0; i < lambda->n; i++) {
+        const double l = psVectorGet(lambda, i);
+
+        if (l > -FLT_EPSILON && l < FLT_EPSILON) {
+            psError(PS_ERR_UNKNOWN, true, "Found almost zero eigenvalue %g.  Returning NULL.\n", l);
+            psFree(lambda); psFree(evectors);
+            psFree(A);
+            psFree(B);
+            return false;
+        }
+    }
+
+    psFree(lambda);
+    psFree(evectors);
+#endif
     if (!psMatrixGJSolve(A, B)) {
         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
@@ -1116 +1155 @@
         result = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64);
         if (!result) {
-            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
         }
         break;
@@ -1253 +1292 @@
         if (!psVectorStats(stats, resid, NULL, mask, maskValue)) {
             psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
-            psFree(resid)
-            psFree(fit)
+            psFree(resid);
+            psFree(fit);
             return false;
         }
@@ -1260 +1299 @@
         double meanValue = psStatsGetValue (stats, meanOption);
         double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        // correct for the number of degrees of freedom absorbed by the fit
+        const int nterm = psPolynomial2DNterm(poly);
+        if (Nkeep <= nterm) {   // we can't estimate the standard deviation
+            psError(PS_ERR_UNKNOWN, true, "Number of data points %d is <= number of degrees of freedom %d",
+                    Nkeep, nterm);
+            psFree(resid);
+            psFree(fit);
+            
+            return false;
+        }
+        stdevValue *= (Nkeep - 1)/(Nkeep - nterm - 1);
 
         psTrace("psLib.math", 5, "Mean is %f\n", meanValue);