The estimation package provides classes to fit some non-linear model to available observations depending on it. These problems are commonly called estimation problems.
The estimation problems considered here are parametric problems where a user-provided model depends on initially unknown scalar parameters and several measurements made on values that depend on the model are available. As examples, one can consider the center and radius of a circle given points approximately lying on a ring, or a satellite orbit given range, range-rate and angular measurements from various ground stations.
One important class of estimation problems is weighted least squares problems. They basically consist in finding the values for some parameters pk such that a cost function J = sum(wiri2) is minimized. The various ri terms represent the deviation ri = mesi - modi between the measurements and the parameterized models. The wi factors are the measurements weights, they are often chosen either all equal to 1.0 or proportional to the inverse of the variance of the measurement type. The solver adjusts the values of the estimated parameters pk which are not bound (i.e. the free parameters). It does not touch the parameters which have been put in a bound state by the user.
The aim of this package is similar to the aim of the optimization package, but the algorithms are entirely different as:
The problem modeling is the most important part for the user. Understanding it is the key to proper use of the package. One interface and two classes are provided for this purpose: EstimationProblem , EstimatedParameter and WeightedMeasurement .
Consider the following example problem: we want to determine the linear trajectory of a sailing ship by performing angular and distance measurements from an observing spot on the shore. The problem model is represented by two equations:
x(t) = x0+(t-t0)vx0
y(t) = y0+(t-t0)vy0
These two equations depend on four parameters (x0, y0, vx0 and vy0). We want to determine these four parameters.
Assuming the observing spot is located at the origin of the coordinates system and that the angular measurements correspond to the angle between the x axis and the line of sight, the theoretical values of the angular measurements at ti and of the distance measurements at tj are modeled as follows:
anglei,theo = atan2(y(ti), x(ti))
distancej,theo = sqrt(x(tj)2+y(tj)2)
The real observations generate a set of measurements values anglei,meas and distancej,meas.
The following class diagram shows one way to solve this problem using the estimation package. The grey elements are already provided by the package whereas the purple elements are developed by the user.
The TrajectoryDeterminationProblem
class holds the linear model
equations x(t) and y(t). It delegate storage of the four parameters x0,
y0, vx0 and vy0 and of the various measurements
anglei,meas and distancej,meas to its base class
SimpleEstimationProblem
. Since the theoretical values of the measurements
anglei,theo and distancej,theo depend on the linear model,
the two classes AngularMeasurement
and DistanceMeasurement
are implemented as internal classes, thus having access to the equations of the
linear model and to the parameters.
Here are the various parts of the TrajectoryDeterminationProblem.java
source file. This example, with an additional main
method is
available here
.
public class TrajectoryDeterminationProblem extends SimpleEstimationProblem { public TrajectoryDeterminationProblem(double t0, double x0Guess, double y0Guess, double vx0Guess, double vy0Guess) { this.t0 = t0; x0 = new EstimatedParameter( "x0", x0Guess); y0 = new EstimatedParameter( "y0", y0Guess); vx0 = new EstimatedParameter("vx0", vx0Guess); vy0 = new EstimatedParameter("vy0", vy0Guess); // inform the base class about the parameters addParameter(x0); addParameter(y0); addParameter(vx0); addParameter(vy0); } public double getX0() { return x0.getEstimate(); } public double getY0() { return y0.getEstimate(); } public double getVx0() { return vx0.getEstimate(); } public double getVy0() { return vy0.getEstimate(); } public void addAngularMeasurement(double wi, double ti, double ai) { // let the base class handle the measurement addMeasurement(new AngularMeasurement(wi, ti, ai)); } public void addDistanceMeasurement(double wi, double ti, double di) { // let the base class handle the measurement addMeasurement(new DistanceMeasurement(wi, ti, di)); } public double x(double t) { return x0.getEstimate() + (t - t0) * vx0.getEstimate(); } public double y(double t) { return y0.getEstimate() + (t - t0) * vy0.getEstimate(); } // measurements internal classes go here private double t0; private EstimatedParameter x0; private EstimatedParameter y0; private EstimatedParameter vx0; private EstimatedParameter vy0; }
serialVersionUID
static fields are present because
the WeightedMeasurement
class implements the
Serializable
interface.
private class AngularMeasurement extends WeightedMeasurement { public AngularMeasurement(double weight, double t, double angle) { super(weight, angle); this.t = t; } public double getTheoreticalValue() { return Math.atan2(y(t), x(t)); } public double getPartial(EstimatedParameter parameter) { double xt = x(t); double yt = y(t); double r = Math.sqrt(xt * xt + yt * yt); double u = yt / (r + xt); double c = 2 * u / (1 + u * u); if (parameter == x0) { return -c; } else if (parameter == vx0) { return -c * t; } else if (parameter == y0) { return c * xt / yt; } else { return c * t * xt / yt; } } private final double t; private static final long serialVersionUID = -5990040582592763282L; }
private class DistanceMeasurement extends WeightedMeasurement { public DistanceMeasurement(double weight, double t, double angle) { super(weight, angle); this.t = t; } public double getTheoreticalValue() { double xt = x(t); double yt = y(t); return Math.sqrt(xt * xt + yt * yt); } public double getPartial(EstimatedParameter parameter) { double xt = x(t); double yt = y(t); double r = Math.sqrt(xt * xt + yt * yt); if (parameter == x0) { return xt / r; } else if (parameter == vx0) { return xt * t / r; } else if (parameter == y0) { return yt / r; } else { return yt * t / r; } } private final double t; private static final long serialVersionUID = 3257286197740459503L; }
estimate
method. Two implementations are already provided by the library:
GaussNewtonEstimator
and
LevenbergMarquardtEstimator
. The first one implements a simple Gauss-Newton
algorithm, which is sufficient when the starting point (initial guess) is close
enough to the solution. The second one implements a more complex Levenberg-Marquardt
algorithm which is more robust when the initial guess is far from the solution.
The following sequence diagram explains roughly what occurs under the hood
in the estimate
method.
Basically, the estimator first retrieves the parameters and the measurements. The estimation loop is based on the gradient of the sum of the squares of the residuals, hence, the estimators get the various partial derivatives of all measurements with respect to all parameters. A new state hopefully globally reducing the residuals is estimated, and the parameters value are updated. This estimation loops stops when either the convergence conditions are met or the maximal number of iterations is exceeded.
One important tuning parameter for weighted least-squares solving is the weight attributed to each measurement. This weights has two purposes:
The weight is a multiplicative factor for the square of the residuals. A common choice is to use the inverse of the variance of the measurements error as the weighting factor for all measurements for one type. On our sailing ship example, we may have a range measurements accuracy of about 1 meter and an angular measurements accuracy of about 0.01 degree, or 1.7 10-4 radians. So we would use w=1.0 for distance measurements weight and w=3 107 for angular measurements weight. If we knew that the measurements quality is bad at tracking start because of measurement system warm-up delay for example, then we would reduce the weight for the first measurements and use for example w=0.1 and w=3 106 respectively, depending on the type.
After a problem has been set up, it is possible to fine tune the
way it will be solved. For example, it may appear the measurements are not
sufficient to get some parameters with sufficient confidence due to observability
problems. It is possible to fix some parameters in order to prevent the solver
from changing them. This is realized by passing true
to the
setBound
method of the parameter.
It is also possible to ignore some measurements by passing true
to the
setIgnored
method of the measurement. A typical use is to