You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/09/25 20:33:34 UTC
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/
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;