MATH-1379: Fix "LoessInterpolator" in case of unevenly spaced samples.

Thanks to Richard Wilkinson.
This commit is contained in:
Gilles Sadowski 2021-07-19 18:40:52 +02:00
parent 30a2593c2c
commit 8968416790
3 changed files with 86 additions and 3 deletions

View File

@ -409,7 +409,8 @@ public class LoessInterpolator
* {@code (right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])}.
* The array will be updated.
*/
private static void updateBandwidthInterval(final double[] xval, final double[] weights,
private static void updateBandwidthInterval(final double[] xval,
final double[] weights,
final int i,
final int[] bandwidthInterval) {
final int left = bandwidthInterval[0];
@ -418,10 +419,13 @@ public class LoessInterpolator
// The right edge should be adjusted if the next point to the right
// is closer to xval[i] than the leftmost point of the current interval
int nextRight = nextNonzero(weights, right);
if (nextRight < xval.length && xval[nextRight] - xval[i] < xval[i] - xval[left]) {
int nextLeft = nextNonzero(weights, bandwidthInterval[0]);
int nextLeft = left;
while (nextRight < xval.length &&
xval[nextRight] - xval[i] < xval[i] - xval[nextLeft]) {
nextLeft = nextNonzero(weights, bandwidthInterval[0]);
bandwidthInterval[0] = nextLeft;
bandwidthInterval[1] = nextRight;
nextRight = nextNonzero(weights, nextRight);
}
}

View File

@ -247,6 +247,82 @@ public class LoessInterpolatorTest {
}
}
// MATH-1379
@Test
public void testFitWithUnevenXSpacing() {
final double[] xval = {
0.1, 0.12, 0.23, 0.4, 0.57,
0.7, 0.87, 1.3, 1.9, 2.2,
2.3, 2.65, 3.0, 3.1, 3.5,
4.6, 4.7, 5.8, 5.95, 6.1
};
final 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
};
// Compare with output from R.
// predict(loess(y ~ x, data.frame(x=xval, y=yval), span=0.35, degree=1, family="symmetric", control=loess.control(iterations=1, surface="direct")))
final double[] yref = {
0.556184894, 0.541907126, 0.455059334, 0.303681477, 0.142126445,
0.002615653, -0.031178445, -0.187124310, -0.405235207, -0.535023851,
-0.706801740, -1.466740294, -2.349248503, -2.596576469, -3.354222419,
0.086206868, 0.320251370, 3.064778450, 3.426179479, 3.783500164
};
final double delta = 1e-8;
// Note R counts all iterations whereas LoessInterpolator robustness iterations are
// in addition to the initial fit.
final LoessInterpolator li = new LoessInterpolator(0.35, 0, 1e-12);
final 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], delta);
}
}
// MATH-1379
@Test
public void testFitWithVaryingWeightsAndUnevenXSpacing() {
final double[] xval = {
0.1, 0.12, 0.23, 0.4, 0.57,
0.7, 0.87, 1.3, 1.9, 2.2,
2.3, 2.65, 3.0, 3.1, 3.5,
4.6, 4.7, 5.8, 5.95, 6.1
};
final 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};
final double[] weights = {
1, 1, 1, 1, 1,
1, 1, 1, 1, 1,
1, 0.8, 0.5, 0.5, 0.8,
0.8, 0.5, 0.5, 0.9, 1
};
// Compare with output from R.
// predict(loess(y ~ x, data.frame(x=xval, y=yval), weights, span=0.35, degree=1, family="symmetric", control=loess.control(iterations=1, surface="direct")))
final double[] yref = {
0.556184894, 0.541907126, 0.455059334, 0.303681477, 0.142126445,
0.002615653, -0.031178445, -0.187124310, -0.406403569, -0.531957113,
-0.669978426, -1.411039850, -2.225022609, -2.463874517, -3.298767758,
-0.506116921, -0.231242991, 2.873755984, 3.295897106, 3.715495321};
final double delta = 1e-8;
// Note R counts all iterations whereas LoessInterpolator robustness iterations are
// in addition to the initial fit.
final LoessInterpolator li = new LoessInterpolator(0.35, 0, 1e-12);
final 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], delta);
}
}
private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
double dx = 2 * AccurateMath.PI / xval.length;
double x = 0;

View File

@ -95,6 +95,9 @@ Caveat:
nightmare was one of the main reasons for creating more focused
components.]
">
<action dev="erans" type="fix" issue="MATH-1379" due-to="Richard Wilkinson">
Fix "LoessInterpolator" (in package "o.a.c.m.legacy.analysis.interpolation").
</action>
<action dev="erans" type="fix" issue="MATH-1617" due-to="Ng Tsz Sum">
Class "BigReal": Fix equality check.
</action>