Browse Source

Merge pull request #685 from diluculo/issue684

Fix incorrect weighted Jacobian in NonlinearObjectiveFunction
pull/697/head
Christoph Ruegg 6 years ago
committed by GitHub
parent
commit
29b22cb2d7
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
  1. 35
      src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs
  2. 6
      src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs

35
src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs

@ -618,5 +618,40 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
}
#endregion Thurber
#region Weighted Nonlinear Regression
// Data from https://www.mathworks.com/help/stats/examples/weighted-nonlinear-regression.html
private Vector<double> PollutionModel(Vector<double> p, Vector<double> x)
{
var y = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
y[i] = p[0] * (1.0 - Math.Exp(-p[1] * x[i]));
}
return y;
}
private Vector<double> PollutionX = new DenseVector(new double[] { 1, 2, 3, 5, 7, 10 });
private Vector<double> PollutionY = new DenseVector(new double[] { 109, 149, 149, 191, 213, 224 });
private Vector<double> PollutionW = new DenseVector(new double[] { 1, 1, 5, 5, 5, 5 });
private Vector<double> PollutionStart = new DenseVector(new double[] { 240, 0.5 });
private Vector<double> PollutionBest = new DenseVector(new double[] { 225.17, 0.40078 });
[Test]
public void PollutionWithWeights()
{
var obj = ObjectiveFunction.NonlinearModel(PollutionModel, PollutionX, PollutionY, PollutionW, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, PollutionStart);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(PollutionBest[i], result.MinimizingPoint[i], 4);
}
}
#endregion Weighted Nonlinear Regression
}
}

6
src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs

@ -329,11 +329,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveFunctions
// if j-th parameter is fixed, set J[i, j] = 0
jacobianValue[i, j] = 0.0;
}
else
else if (Weights != null)
{
jacobianValue[i, j] = (Weights == null)
? jacobianValue[i, j]
: jacobianValue[i, j] * L[j];
jacobianValue[i, j] = jacobianValue[i, j] * L[i];
}
}
}

Loading…
Cancel
Save