mirror of
https://github.com/apache/commons-math.git
synced 2025-03-08 01:20:01 +00:00
fixed wrong results in Loess interpolator
also added a way to set weights for smoothing JIRA: MATH-296 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@818942 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
97fb25eb35
commit
d98e430b04
src
main/java/org/apache/commons/math
site/xdoc
test/java/org/apache/commons/math/analysis/interpolation
@ -188,6 +188,8 @@ public class MessagesResources_fr
|
||||
"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, " +
|
||||
|
||||
@ -47,6 +47,9 @@ public class LoessInterpolator
|
||||
/** 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;
|
||||
|
||||
@ -69,21 +72,34 @@ public class LoessInterpolator
|
||||
*/
|
||||
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 @@ public class LoessInterpolator
|
||||
* {@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 @@ public class LoessInterpolator
|
||||
robustnessIters);
|
||||
}
|
||||
this.robustnessIters = robustnessIters;
|
||||
this.accuracy = accuracy;
|
||||
}
|
||||
|
||||
/**
|
||||
@ -135,10 +178,11 @@ public class LoessInterpolator
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 @@ public class LoessInterpolator
|
||||
* <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 @@ public class LoessInterpolator
|
||||
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 @@ public class LoessInterpolator
|
||||
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 @@ public class LoessInterpolator
|
||||
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,19 +325,48 @@ public class LoessInterpolator
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
/**
|
||||
* 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
|
||||
@ -336,18 +411,14 @@ public class LoessInterpolator
|
||||
* 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);
|
||||
}
|
||||
}
|
||||
|
||||
@ -44,6 +44,10 @@ The <action> type attribute can be add,update,fix,remove.
|
||||
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>
|
||||
|
||||
@ -219,6 +219,54 @@ public class LoessInterpolatorTest {
|
||||
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;
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user