Return-Path: Delivered-To: apmail-commons-commits-archive@locus.apache.org Received: (qmail 71472 invoked from network); 5 Feb 2008 18:09:39 -0000 Received: from hermes.apache.org (HELO mail.apache.org) (140.211.11.2) by minotaur.apache.org with SMTP; 5 Feb 2008 18:09:39 -0000 Received: (qmail 25248 invoked by uid 500); 5 Feb 2008 18:09:30 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 25179 invoked by uid 500); 5 Feb 2008 18:09:29 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 25170 invoked by uid 99); 5 Feb 2008 18:09:29 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Tue, 05 Feb 2008 10:09:29 -0800 X-ASF-Spam-Status: No, hits=-100.0 required=10.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.3] (HELO eris.apache.org) (140.211.11.3) by apache.org (qpsmtpd/0.29) with ESMTP; Tue, 05 Feb 2008 18:09:08 +0000 Received: by eris.apache.org (Postfix, from userid 65534) id 537251A9832; Tue, 5 Feb 2008 10:09:16 -0800 (PST) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r618726 - in /commons/proper/math/trunk: src/site/resources/userguide/ xdocs/userguide/ Date: Tue, 05 Feb 2008 18:09:15 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.8 Message-Id: <20080205180916.537251A9832@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org Author: luc Date: Tue Feb 5 10:09:06 2008 New Revision: 618726 URL: http://svn.apache.org/viewvc?rev=618726&view=rev Log: completely rewrote estimation package documentation with downloadable example and explanation diagrams Added: commons/proper/math/trunk/src/site/resources/userguide/ commons/proper/math/trunk/src/site/resources/userguide/TrajectoryDeterminationProblem.java (with props) commons/proper/math/trunk/src/site/resources/userguide/estimation-class-diagram.png (with props) commons/proper/math/trunk/src/site/resources/userguide/estimation-sequence-diagram.png (with props) Modified: commons/proper/math/trunk/xdocs/userguide/estimation.xml Added: commons/proper/math/trunk/src/site/resources/userguide/TrajectoryDeterminationProblem.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/resources/userguide/TrajectoryDeterminationProblem.java?rev=618726&view=auto ============================================================================== --- commons/proper/math/trunk/src/site/resources/userguide/TrajectoryDeterminationProblem.java (added) +++ commons/proper/math/trunk/src/site/resources/userguide/TrajectoryDeterminationProblem.java Tue Feb 5 10:09:06 2008 @@ -0,0 +1,187 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +import org.apache.commons.math.estimation.EstimationException; +import org.apache.commons.math.estimation.EstimatedParameter; +import org.apache.commons.math.estimation.EstimationProblem; +import org.apache.commons.math.estimation.LevenbergMarquardtEstimator; +import org.apache.commons.math.estimation.SimpleEstimationProblem; +import org.apache.commons.math.estimation.WeightedMeasurement; + +public class TrajectoryDeterminationProblem extends SimpleEstimationProblem { + + public static void main(String[] args) { + try { + TrajectoryDeterminationProblem problem = + new TrajectoryDeterminationProblem(0.0, 100.0, 800.0, 1.0, 0.0); + + double[][] distances = { + { 0.0, 806.5849 }, { 20.0, 796.8148 }, { 40.0, 791.0833 }, { 60.0, 789.6712 }, + { 80.0, 793.1334 }, { 100.0, 797.7248 }, { 120.0, 803.2785 }, { 140.0, 813.4939 }, + { 160.0, 826.9295 }, { 180.0, 844.0640 }, { 200.0, 863.3829 }, { 220.0, 883.3143 }, + { 240.0, 908.6867 }, { 260.0, 934.8561 }, { 280.0, 964.0730 }, { 300.0, 992.1033 }, + { 320.0, 1023.998 }, { 340.0, 1057.439 }, { 360.0, 1091.912 }, { 380.0, 1125.968 }, + { 400.0, 1162.789 }, { 420.0, 1201.517 }, { 440.0, 1239.176 }, { 460.0, 1279.347 } }; + for (int i = 0; i < distances.length; ++i) { + problem.addDistanceMeasurement(1.0, distances[i][0], distances[i][1]); + }; + + double[][] angles = { + { 10.0, 1.415423 }, { 30.0, 1.352643 }, { 50.0, 1.289290 }, { 70.0, 1.225249 }, + { 90.0, 1.161203 }, {110.0, 1.098538 }, {130.0, 1.036263 }, {150.0, 0.976052 }, + {170.0, 0.917921 }, {190.0, 0.861830 }, {210.0, 0.808237 }, {230.0, 0.757043 }, + {250.0, 0.708650 }, {270.0, 0.662949 }, {290.0, 0.619903 }, {310.0, 0.579160 }, + {330.0, 0.541033 }, {350.0, 0.505590 }, {370.0, 0.471746 }, {390.0, 0.440155 }, + {410.0, 0.410522 }, {430.0, 0.382701 }, {450.0, 0.356957 }, {470.0, 0.332400 } }; + for (int i = 0; i < angles.length; ++i) { + problem.addAngularMeasurement(3.0e7, angles[i][0], angles[i][1]); + }; + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + System.out.println("initial position: " + problem.getX0() + " " + problem.getY0()); + System.out.println("velocity: " + problem.getVx0() + " " + problem.getVy0()); + + } catch (EstimationException ee) { + System.err.println(ee.getMessage()); + } + } + + 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(); + } + + 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; + + } + + private double t0; + private EstimatedParameter x0; + private EstimatedParameter y0; + private EstimatedParameter vx0; + private EstimatedParameter vy0; + +} Propchange: commons/proper/math/trunk/src/site/resources/userguide/TrajectoryDeterminationProblem.java ------------------------------------------------------------------------------ svn:eol-style = native Added: commons/proper/math/trunk/src/site/resources/userguide/estimation-class-diagram.png URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/resources/userguide/estimation-class-diagram.png?rev=618726&view=auto ============================================================================== Binary file - no diff available. Propchange: commons/proper/math/trunk/src/site/resources/userguide/estimation-class-diagram.png ------------------------------------------------------------------------------ svn:mime-type = image/png Added: commons/proper/math/trunk/src/site/resources/userguide/estimation-sequence-diagram.png URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/resources/userguide/estimation-sequence-diagram.png?rev=618726&view=auto ============================================================================== Binary file - no diff available. Propchange: commons/proper/math/trunk/src/site/resources/userguide/estimation-sequence-diagram.png ------------------------------------------------------------------------------ svn:mime-type = image/png Modified: commons/proper/math/trunk/xdocs/userguide/estimation.xml URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/xdocs/userguide/estimation.xml?rev=618726&r1=618725&r2=618726&view=diff ============================================================================== --- commons/proper/math/trunk/xdocs/userguide/estimation.xml (original) +++ commons/proper/math/trunk/xdocs/userguide/estimation.xml Tue Feb 5 10:09:06 2008 @@ -1,22 +1,22 @@ + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +--> + @@ -29,86 +29,336 @@

- 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(wi ri2) - 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 estimation package provides classes to fit some non-linear + model to available observations depending on it. These + problems are commonly called estimation problems.

- The aim of this package is similar to the aim of the optimization package, but the - algorithms are entirely differents as: + 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:

  • - they need the partial derivatives of the measurements - with respect to the free parameters + they need the partial derivatives of the measurements with + respect to the free parameters
  • - they are residuals based instead of generic cost functions based + they are residuals based instead of generic cost functions + based

- - - - The - org.apache.commons.math.estimation.EstimationProblem interface is a very - simple container packing together parameters and measurements. - + +

- The - org.apache.commons.math.estimation.EstimatedParameter class to represent each - estimated parameter. The parameters are set up with a guessed value which will be - adjusted by the solver until a best fit is achieved. It is possible to change which - parameters are modified and which are preserved thanks to a bound property. Such - settings are often needed by expert users to handle contingency cases with very - low observability. + 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.

- -

- The user extends the - org.apache.commons.math.estimation.WeightedMeasurement abstract class to define its - own measurements. Each measurement types should have its own implementing class, for - example in the satellite example above , the user should define three classes, one - for range measurements, one for range-rates measurements and one for angular measurements. - Each measurement would correspond to an instance of the appropriate class, set up with - the date, a reference to the ground station, the weight and the measured value. -

- - + 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: +

- The package provides two common - org.apache.commons.math.estimation.Estimator implementations to solve weighted - least squares problems. The first one is based on the - Gauss-Newton method. - The second one is based on the - Levenberg-Marquardt - method. The first one is best suited when a good approximation of the parameters is known while the second one - is more robust and can handle starting points far from the solution. + 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. +

+
First, the general setup of the class: declarations, fields, constructor, setters and getters: + +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; + +} + +
+
The two specialized measurements class are simple internal classes that + implement the equation for their respective measurement type, using the + enclosing class to get the parameters references and the linear models x(t) + and y(t). The 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; + + } + +
+ + +

+ Solving the problem is simply a matter of choosing an implementation + of the + Estimator interface and to pass the problem instance to its 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: +

+
    +
  • fixing unit problems when combining different types of measurements
    • +
  • adjusting the influence of good or bad measurements on the solution
    • +
+

+ 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 +

    +
  1. + perform a first determination with all parameters, to check each measurement + residual after convergence (i.e. to compute the difference between the + measurement and its theoretical value as computed from the estimated parameters), +
    2. +
  2. + compute standard deviation for the measurements samples (one sample for each + measurements type) +
    3. +
  3. + ignore measurements whose residual are above some threshold (for example three + time the standard deviation on the residuals) assuming they correspond to + bad measurements, +
    4. +
  4. + perform another determination on the reduced measurements set. +
    5. +
+

+