commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r818942 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ main/java/org/apache/commons/math/analysis/interpolation/ site/xdoc/ test/java/org/apache/commons/math/analysis/interpolation/
Date Fri, 25 Sep 2009 18:33:34 GMT
Author: luc
Date: Fri Sep 25 18:33:34 2009
New Revision: 818942

URL: http://svn.apache.org/viewvc?rev=818942&view=rev
Log:
fixed wrong results in Loess interpolator
also added a way to set weights for smoothing
JIRA: MATH-296

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/MessagesResources_fr.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/MessagesResources_fr.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/MessagesResources_fr.java?rev=818942&r1=818941&r2=818942&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/MessagesResources_fr.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/MessagesResources_fr.java
Fri Sep 25 18:33:34 2009
@@ -188,6 +188,8 @@
       "toutes les abscisses doivent \u00eatre des nombres r\u00e9els finis, mais l''abscisse
{0} vaut {1}" },
     { "all ordinatae must be finite real numbers, but {0}-th is {1}",
       "toutes les ordonn\u00e9es doivent \u00eatre des nombres r\u00e9els finis, mais l''ordonn\u00e9e
{0} vaut {1}" },
+    { "all weights must be finite real numbers, but {0}-th is {1}",
+      "tous les poids doivent \u00eatre des nombres r\u00e9els finis, mais le poids {0} vaut
{1}" },
     { "the abscissae array must be sorted in a strictly increasing order, " +
       "but the {0}-th element is {1} whereas {2}-th is {3}",
       "les abscisses doivent \u00eatre en ordre strictement croissant, " +

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java?rev=818942&r1=818941&r2=818942&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
Fri Sep 25 18:33:34 2009
@@ -47,6 +47,9 @@
     /** Default value of the number of robustness iterations. */
     public static final int DEFAULT_ROBUSTNESS_ITERS = 2;
 
+    /** Default value for accuracy. */
+    public static final double DEFAULT_ACCURACY = 1e-12;
+
     /** serializable version identifier. */
     private static final long serialVersionUID = 5204927143605193821L;
 
@@ -70,20 +73,33 @@
     private final int robustnessIters;
 
     /**
+     * If the median residual at a certain robustness iteration
+     * is less than this amount, no more iterations are done.
+     */
+    private final double accuracy;
+
+    /**
      * Constructs a new {@link LoessInterpolator}
-     * with a bandwidth of {@link #DEFAULT_BANDWIDTH} and
-     * {@link #DEFAULT_ROBUSTNESS_ITERS} robustness iterations.
-     * See {@link #LoessInterpolator(double, int)} for an explanation of
+     * with a bandwidth of {@link #DEFAULT_BANDWIDTH},
+     * {@link #DEFAULT_ROBUSTNESS_ITERS} robustness iterations
+     * and an accuracy of {#link #DEFAULT_ACCURACY}.
+     * See {@link #LoessInterpolator(double, int, double)} for an explanation of
      * the parameters.
      */
     public LoessInterpolator() {
         this.bandwidth = DEFAULT_BANDWIDTH;
         this.robustnessIters = DEFAULT_ROBUSTNESS_ITERS;
+        this.accuracy = DEFAULT_ACCURACY;
     }
 
     /**
      * Constructs a new {@link LoessInterpolator}
      * with given bandwidth and number of robustness iterations.
+     * <p>
+     * Calling this constructor is equivalent to calling {link {@link
+     * #LoessInterpolator(double, int, double) LoessInterpolator(bandwidth,
+     * robustnessIters, LoessInterpolator.DEFAULT_ACCURACY)}
+     * </p>
      *
      * @param bandwidth  when computing the loess fit at
      * a particular point, this fraction of source points closest
@@ -97,8 +113,34 @@
      * {@link #DEFAULT_ROBUSTNESS_ITERS}.
      * @throws MathException if bandwidth does not lie in the interval [0,1]
      * or if robustnessIters is negative.
+     * @see #LoessInterpolator(double, int, double)
      */
     public LoessInterpolator(double bandwidth, int robustnessIters) throws MathException
{
+        this(bandwidth, robustnessIters, DEFAULT_ACCURACY);
+    }
+
+    /**
+     * Constructs a new {@link LoessInterpolator}
+     * with given bandwidth, number of robustness iterations and accuracy.
+     *
+     * @param bandwidth  when computing the loess fit at
+     * a particular point, this fraction of source points closest
+     * to the current point is taken into account for computing
+     * a least-squares regression.</br>
+     * A sensible value is usually 0.25 to 0.5, the default value is
+     * {@link #DEFAULT_BANDWIDTH}.
+     * @param robustnessIters This many robustness iterations are done.</br>
+     * A sensible value is usually 0 (just the initial fit without any
+     * robustness iterations) to 4, the default value is
+     * {@link #DEFAULT_ROBUSTNESS_ITERS}.
+     * @param accuracy If the median residual at a certain robustness iteration
+     * is less than this amount, no more iterations are done.
+     * @throws MathException if bandwidth does not lie in the interval [0,1]
+     * or if robustnessIters is negative.
+     * @see #LoessInterpolator(double, int)
+     * @since 2.1
+     */
+    public LoessInterpolator(double bandwidth, int robustnessIters, double accuracy) throws
MathException {
         if (bandwidth < 0 || bandwidth > 1) {
             throw new MathException("bandwidth must be in the interval [0,1], but got {0}",
                                     bandwidth);
@@ -110,6 +152,7 @@
                                     robustnessIters);
         }
         this.robustnessIters = robustnessIters;
+        this.accuracy = accuracy;
     }
 
     /**
@@ -135,10 +178,11 @@
     }
 
     /**
-     * Compute a loess fit on the data at the original abscissae.
+     * Compute a weighted loess fit on the data at the original abscissae.
      *
      * @param xval the arguments for the interpolation points
      * @param yval the values for the interpolation points
+     * @param weights point weights: coefficients by which the robustness weight of a point
is multiplied
      * @return values of the loess fit at corresponding original abscissae
      * @throws MathException if some of the following conditions are false:
      * <ul>
@@ -146,8 +190,9 @@
      * <li> The arguments are in a strictly increasing order</li>
      * <li> All arguments and values are finite real numbers</li>
      * </ul>
+     * @since 2.1
      */
-    public final double[] smooth(final double[] xval, final double[] yval)
+    public final double[] smooth(final double[] xval, final double[] yval, final double[]
weights)
             throws MathException {
         if (xval.length != yval.length) {
             throw new MathException(
@@ -163,8 +208,9 @@
             throw new MathException("Loess expects at least 1 point");
         }
 
-        checkAllFiniteReal(xval, true);
-        checkAllFiniteReal(yval, false);
+        checkAllFiniteReal(xval, "all abscissae must be finite real numbers, but {0}-th is
{1}");
+        checkAllFiniteReal(yval, "all ordinatae must be finite real numbers, but {0}-th is
{1}");
+        checkAllFiniteReal(weights, "all weights must be finite real numbers, but {0}-th
is {1}");
 
         checkStrictlyIncreasing(xval);
 
@@ -237,7 +283,7 @@
                     final double xk   = xval[k];
                     final double yk   = yval[k];
                     final double dist = (k < i) ? x - xk : xk - x;
-                    final double w    = tricube(dist * denom) * robustnessWeights[k];
+                    final double w    = tricube(dist * denom) * robustnessWeights[k] * weights[k];
                     final double xkw  = xk * w;
                     sumWeights += w;
                     sumX += xkw;
@@ -252,7 +298,7 @@
                 final double meanXSquared = sumXSquared / sumWeights;
 
                 final double beta;
-                if (meanXSquared == meanX * meanX) {
+                if (Math.sqrt(Math.abs(meanXSquared - meanX * meanX)) < accuracy) {
                     beta = 0;
                 } else {
                     beta = (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX);
@@ -279,13 +325,18 @@
             Arrays.sort(sortedResiduals);
             final double medianResidual = sortedResiduals[n / 2];
 
-            if (medianResidual == 0) {
+            if (Math.abs(medianResidual) < accuracy) {
                 break;
             }
 
             for (int i = 0; i < n; ++i) {
                 final double arg = residuals[i] / (6 * medianResidual);
-                robustnessWeights[i] = (arg >= 1) ? 0 : Math.pow(1 - arg * arg, 2);
+                if (arg >= 1) {
+                    robustnessWeights[i] = 0;
+                } else {
+                    final double w = 1 - arg * arg;
+                    robustnessWeights[i] = w * w;
+                }
             }
         }
 
@@ -293,6 +344,30 @@
     }
 
     /**
+     * Compute a loess fit on the data at the original abscissae.
+     *
+     * @param xval the arguments for the interpolation points
+     * @param yval the values for the interpolation points
+     * @return values of the loess fit at corresponding original abscissae
+     * @throws MathException if some of the following conditions are false:
+     * <ul>
+     * <li> Arguments and values are of the same size that is greater than zero</li>
+     * <li> The arguments are in a strictly increasing order</li>
+     * <li> All arguments and values are finite real numbers</li>
+     * </ul>
+     */
+    public final double[] smooth(final double[] xval, final double[] yval)
+            throws MathException {
+
+        final double[] unitWeights = new double[xval.length];
+        Arrays.fill(unitWeights, 1.0);
+
+        return smooth(xval, yval, unitWeights);
+
+    }
+
+
+    /**
      * Given an index interval into xval that embraces a certain number of
      * points closest to xval[i-1], update the interval so that it embraces
      * the same number of points closest to xval[i]
@@ -336,18 +411,14 @@
      * Check that all elements of an array are finite real numbers.
      *
      * @param values the values array
-     * @param isAbscissae if true, elements are abscissae otherwise they are ordinatae
-     * @throws MathException if one of the values is not
-     *         a finite real number
+     * @param pattern pattern of the error message
+     * @throws MathException if one of the values is not a finite real number
      */
-    private static void checkAllFiniteReal(final double[] values, final boolean isAbscissae)
+    private static void checkAllFiniteReal(final double[] values, final String pattern)
         throws MathException {
         for (int i = 0; i < values.length; i++) {
             final double x = values[i];
             if (Double.isInfinite(x) || Double.isNaN(x)) {
-                final String pattern = isAbscissae ?
-                        "all abscissae must be finite real numbers, but {0}-th is {1}" :
-                        "all ordinatae must be finite real numbers, but {0}-th is {1}";
                 throw new MathException(pattern, i, x);
             }
         }

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=818942&r1=818941&r2=818942&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Sep 25 18:33:34 2009
@@ -44,6 +44,10 @@
         interface contract.  Added getGeneratorUpperBounds method to
         EmpiricalDistributionImpl providing previous behavior.
       </action>
+      <action dev="luc" tyoe="fix" issue="MATH-296" due-to="Eugene Kirpichov">
+        Fixed wrong results on Loess interpolation, also added a way to set weights
+        for smoothing
+      </action>
       <action dev="luc" type="fix" issue="MATH-293" due-to="Benjamin McCann">
         Fixed a OutOfBoundException in simplex solver when some constraints are tight.
       </action>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java?rev=818942&r1=818941&r2=818942&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java
Fri Sep 25 18:33:34 2009
@@ -219,6 +219,54 @@
         new LoessInterpolator(1.1, 3);
     }
 
+    @Test
+    public void testMath296withoutWeights() throws MathException {
+        double[] xval = {
+                0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
+                 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0};
+        double[] yval = {
+                0.47, 0.48, 0.55, 0.56, -0.08, -0.04, -0.07, -0.07,
+                -0.56, -0.46, -0.56, -0.52, -3.03, -3.08, -3.09,
+                -3.04, 3.54, 3.46, 3.36, 3.35};
+        // Output from R, rounded to .001
+        double[] yref = {
+                0.461, 0.499, 0.541, 0.308, 0.175, -0.042, -0.072,
+                -0.196, -0.311, -0.446, -0.557, -1.497, -2.133,
+                -3.08, -3.09, -0.621, 0.982, 3.449, 3.389, 3.336
+        };
+        LoessInterpolator li = new LoessInterpolator(0.3, 4, 1e-12);
+        double[] res = li.smooth(xval, yval);
+        Assert.assertEquals(xval.length, res.length);
+        for(int i = 0; i < res.length; ++i) {
+            Assert.assertEquals(yref[i], res[i], 0.02);
+        }
+    }
+
+    @Test
+    public void testMath296withWeights() throws MathException {
+        double[] xval = {
+                0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
+                 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0};
+        double[] yval = {
+                0.47, 0.48, 0.55, 0.56, -0.08, -0.04, -0.07, -0.07,
+                -0.56, -0.46, -0.56, -0.52, -3.03, -3.08, -3.09,
+                -3.04, 3.54, 3.46, 3.36, 3.35};
+        double[] weights = {
+                1,1,1,1,1,1,1,1,1,1,
+                1,1,0,0,1,1,0,0,1,1};
+        // Output from R, rounded to .001
+        double[] yref = {
+                0.478, 0.492, 0.484, 0.320, 0.179, -0.003, -0.088, -0.209,
+                -0.327, -0.455, -0.518, -0.537, -1.492, -2.115, -3.09, -3.04,
+                -3.0, 0.155, 1.752, 3.35};
+        LoessInterpolator li = new LoessInterpolator(0.3, 4, 1e-12);
+        double[] res = li.smooth(xval, yval,weights);
+        Assert.assertEquals(xval.length, res.length);
+        for(int i = 0; i < res.length; ++i) {
+            Assert.assertEquals(yref[i], res[i], 0.05);
+        }
+    }
+
     private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise)
{
         double dx = 2 * Math.PI / xval.length;
         double x = 0;



Mime
View raw message