You are viewing a plain text version of this content. The canonical link for it is here.
Posted to issues@commons.apache.org by "Mick (JIRA)" <ji...@apache.org> on 2008/03/20 23:57:24 UTC

[jira] Created: (MATH-199) exception in LevenbergMarquardtEstimator

exception in LevenbergMarquardtEstimator
----------------------------------------

                 Key: MATH-199
                 URL: https://issues.apache.org/jira/browse/MATH-199
             Project: Commons Math
          Issue Type: Bug
    Affects Versions: 1.2
         Environment: Windows XP
Java 6

            Reporter: Mick


I get this exception:

Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
       at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
       at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
       at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
       at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
on the code below.

The exception does not occur all the weights in the quadraticFitter are 0.0;


---------------------------------------------------------------------------------------------

package quadraticFitterProblem;

import org.apache.commons.math.estimation.EstimationException;
import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
//import org.apache.commons.math.estimation.WeightedMeasurement;

import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;

public class QuadraticFitterProblem {

       private QuadraticFitter quadraticFitter;

       public QuadraticFitterProblem() {
         // create the uninitialized fitting problem
         quadraticFitter = new QuadraticFitter();

         quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
         quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
         quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
         quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
         quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);

         try {
           // solve the problem, using a Levenberg-Marquardt algorithm with
default settings
           LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
           //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
           estimator.estimate(quadraticFitter);

         } catch (EstimationException ee) {
               System.err.println(ee.getMessage());
         }
       }

       /**
        * @param args
        *
        */
       public static void main(String[] args) {

                       new QuadraticFitterProblem();
                       System.out.println ("Done.");
       }

}

----------------------------------------------------------------------------------------------
import org.apache.commons.math.estimation.EstimatedParameter;
//import org.apache.commons.math.estimation.EstimationException;
//import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
import org.apache.commons.math.estimation.SimpleEstimationProblem;
import org.apache.commons.math.estimation.WeightedMeasurement;

public class QuadraticFitter extends SimpleEstimationProblem {

       // y = a x<sup>2</sup> + b x + c
   private EstimatedParameter a;
   private EstimatedParameter b;
   private EstimatedParameter c;

   /**
    * constructor
    *
    *Fitter for a quadratic model to a sample of 2D points.
    * <p>The model is y(x) = a x<sup>2</sup> + b x + c
    * its three parameters of the model are a, b and c.</p>
    */
   public QuadraticFitter() {

       // three parameters of the model
       a = new EstimatedParameter("a", 0.0);
       b = new EstimatedParameter("b", 0.0);
       c = new EstimatedParameter("c", 0.0);

       // provide the parameters to the base class which
       // implements the getAllParameters and getUnboundParameters methods
       addParameter(a);
       addParameter(b);
       addParameter(c);
   }

   /**
    * Add a sample point
    *
    * @param x abscissa
    * @param y ordinate
    * @param w weight
    */
   public void addPoint(double x, double y, double w) {
       addMeasurement(new LocalMeasurement(x, y, w));
   }

   /**
    * Get the value of the quadratic coefficient.
    *
    * @return the value of a for the quadratic model
    * y = a x<sup>2</sup> + b x + c
    */
   public double getA() {
       return a.getEstimate();
   }

   /**
    * Get the value of the linear coefficient.
    *
    * @return the value of b for the quadratic model
    * y = a x<sup>2</sup> + b x + c
    */
   public double getB() {
       return b.getEstimate();
   }

   /**
    * Get the value of the constant coefficient.
    *
    * @return the value of ac for the quadratic model
    * y = a x<sup>2</sup> + b x + c
    */
   public double getC() {
       return c.getEstimate();
   }

   /**
    * Get the theoretical value of the model for some x.
    * <p>The theoretical value is the value computed using
    * the current state of the problem parameters.</p>
    *
    * Note the use of Hörner's method (synthetic division) for
evaluating polynomials,
    * (more efficient)
    *
    * @param x explanatory variable
    * @return the theoretical value y = a x<sup>2</sup> + b x + c
    */
   public double theoreticalValue(double x) {
       //System.out.println ("x = " + x + "  a.getEstimate() = " +
a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
c.getEstimate() = " + c.getEstimate());
       return ( (a.getEstimate() * x + b.getEstimate() ) * x +
c.getEstimate());
   }

   /**
    * Get the partial derivative of the theoretical value
    * of the model for some x.
    * <p>The derivative is computed using
    * the current state of the problem parameters.</p>
    *
    * @param x explanatory variable
    * @param parameter estimated parameter (either a, b, or c)
    * @return the partial derivative dy/dp
    */
   private double partial(double x, EstimatedParameter parameter) {
       // since we know the only parameters are a, b and c in this
       // class we simply use "==" for efficiency
       if (parameter == a) {
           return x * x;
       } else if (parameter == b) {
           return x;
       } else {
           return 1.0;
       }

   }


   /** Internal measurements class.
    * <p>The measurement is the y value for a fixed specified x.</p>
    */
   private class LocalMeasurement extends WeightedMeasurement {

       static final long serialVersionUID = 1;

       private final double x;

       // constructor
       public LocalMeasurement(double x, double y, double w) {
           super(w, y);
           this.x = x;
       }

       public double getTheoreticalValue() {
           // the value is provided by the model for the local x
           return theoreticalValue(x);
       }

       public double getPartial(EstimatedParameter parameter) {
           // the value is provided by the model for the local x
           return partial(x, parameter);
       }

   }

 }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.


[jira] Commented: (MATH-199) exception in LevenbergMarquardtEstimator

Posted by "Mick (JIRA)" <ji...@apache.org>.
    [ https://issues.apache.org/jira/browse/MATH-199?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12581577#action_12581577 ] 

Mick commented on MATH-199:
---------------------------

That being the case, could the commons math library define a constant such that its inverse it  equal to Double.MAX_VALUE?
I could then use this as a limit.


> exception in LevenbergMarquardtEstimator
> ----------------------------------------
>
>                 Key: MATH-199
>                 URL: https://issues.apache.org/jira/browse/MATH-199
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 1.2
>         Environment: Windows XP
> Java 6
>            Reporter: Mick
>
> I get this exception:
> Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
>        at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
>        at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
> on the code below.
> The exception does not occur all the weights in the quadraticFitter are 0.0;
> ---------------------------------------------------------------------------------------------
> package quadraticFitterProblem;
> import org.apache.commons.math.estimation.EstimationException;
> import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> //import org.apache.commons.math.estimation.WeightedMeasurement;
> import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;
> public class QuadraticFitterProblem {
>        private QuadraticFitter quadraticFitter;
>        public QuadraticFitterProblem() {
>          // create the uninitialized fitting problem
>          quadraticFitter = new QuadraticFitter();
>          quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
>          quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
>          quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
>          quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
>          quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
>          try {
>            // solve the problem, using a Levenberg-Marquardt algorithm with
> default settings
>            LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
>            //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
>            estimator.estimate(quadraticFitter);
>          } catch (EstimationException ee) {
>                System.err.println(ee.getMessage());
>          }
>        }
>        /**
>         * @param args
>         *
>         */
>        public static void main(String[] args) {
>                        new QuadraticFitterProblem();
>                        System.out.println ("Done.");
>        }
> }
> ----------------------------------------------------------------------------------------------
> import org.apache.commons.math.estimation.EstimatedParameter;
> //import org.apache.commons.math.estimation.EstimationException;
> //import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> import org.apache.commons.math.estimation.SimpleEstimationProblem;
> import org.apache.commons.math.estimation.WeightedMeasurement;
> public class QuadraticFitter extends SimpleEstimationProblem {
>        // y = a x<sup>2</sup> + b x + c
>    private EstimatedParameter a;
>    private EstimatedParameter b;
>    private EstimatedParameter c;
>    /**
>     * constructor
>     *
>     *Fitter for a quadratic model to a sample of 2D points.
>     * <p>The model is y(x) = a x<sup>2</sup> + b x + c
>     * its three parameters of the model are a, b and c.</p>
>     */
>    public QuadraticFitter() {
>        // three parameters of the model
>        a = new EstimatedParameter("a", 0.0);
>        b = new EstimatedParameter("b", 0.0);
>        c = new EstimatedParameter("c", 0.0);
>        // provide the parameters to the base class which
>        // implements the getAllParameters and getUnboundParameters methods
>        addParameter(a);
>        addParameter(b);
>        addParameter(c);
>    }
>    /**
>     * Add a sample point
>     *
>     * @param x abscissa
>     * @param y ordinate
>     * @param w weight
>     */
>    public void addPoint(double x, double y, double w) {
>        addMeasurement(new LocalMeasurement(x, y, w));
>    }
>    /**
>     * Get the value of the quadratic coefficient.
>     *
>     * @return the value of a for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getA() {
>        return a.getEstimate();
>    }
>    /**
>     * Get the value of the linear coefficient.
>     *
>     * @return the value of b for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getB() {
>        return b.getEstimate();
>    }
>    /**
>     * Get the value of the constant coefficient.
>     *
>     * @return the value of ac for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getC() {
>        return c.getEstimate();
>    }
>    /**
>     * Get the theoretical value of the model for some x.
>     * <p>The theoretical value is the value computed using
>     * the current state of the problem parameters.</p>
>     *
>     * Note the use of Hörner's method (synthetic division) for
> evaluating polynomials,
>     * (more efficient)
>     *
>     * @param x explanatory variable
>     * @return the theoretical value y = a x<sup>2</sup> + b x + c
>     */
>    public double theoreticalValue(double x) {
>        //System.out.println ("x = " + x + "  a.getEstimate() = " +
> a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
> c.getEstimate() = " + c.getEstimate());
>        return ( (a.getEstimate() * x + b.getEstimate() ) * x +
> c.getEstimate());
>    }
>    /**
>     * Get the partial derivative of the theoretical value
>     * of the model for some x.
>     * <p>The derivative is computed using
>     * the current state of the problem parameters.</p>
>     *
>     * @param x explanatory variable
>     * @param parameter estimated parameter (either a, b, or c)
>     * @return the partial derivative dy/dp
>     */
>    private double partial(double x, EstimatedParameter parameter) {
>        // since we know the only parameters are a, b and c in this
>        // class we simply use "==" for efficiency
>        if (parameter == a) {
>            return x * x;
>        } else if (parameter == b) {
>            return x;
>        } else {
>            return 1.0;
>        }
>    }
>    /** Internal measurements class.
>     * <p>The measurement is the y value for a fixed specified x.</p>
>     */
>    private class LocalMeasurement extends WeightedMeasurement {
>        static final long serialVersionUID = 1;
>        private final double x;
>        // constructor
>        public LocalMeasurement(double x, double y, double w) {
>            super(w, y);
>            this.x = x;
>        }
>        public double getTheoreticalValue() {
>            // the value is provided by the model for the local x
>            return theoreticalValue(x);
>        }
>        public double getPartial(EstimatedParameter parameter) {
>            // the value is provided by the model for the local x
>            return partial(x, parameter);
>        }
>    }
>  }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.


[jira] Commented: (MATH-199) exception in LevenbergMarquardtEstimator

Posted by "Mick (JIRA)" <ji...@apache.org>.
    [ https://issues.apache.org/jira/browse/MATH-199?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12581536#action_12581536 ] 

Mick commented on MATH-199:
---------------------------

Is it accurate to say that the problem is caused by the 4.4E-323 in the error term?
(ie, the 4.4E-323 causes one of the transformed columns to have values of Infinity or NaN?

If I eliminate those error terms, e.g.:
	  quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
	  quadraticFitter.addPoint (1,  -2.5581184967730577, 0.5);
	  quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
	  quadraticFitter.addPoint (3,  -1.9122489313410047, 0.5);
	  quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
then, indeed, the error does not occur.

I have a concern that, while it is true that  4.4E-323 is a number of extremely small magnitude, it is a vaild double, but using it causes an exception.
:Perhaps I am not seeing this correctly.



> exception in LevenbergMarquardtEstimator
> ----------------------------------------
>
>                 Key: MATH-199
>                 URL: https://issues.apache.org/jira/browse/MATH-199
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 1.2
>         Environment: Windows XP
> Java 6
>            Reporter: Mick
>
> I get this exception:
> Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
>        at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
>        at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
> on the code below.
> The exception does not occur all the weights in the quadraticFitter are 0.0;
> ---------------------------------------------------------------------------------------------
> package quadraticFitterProblem;
> import org.apache.commons.math.estimation.EstimationException;
> import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> //import org.apache.commons.math.estimation.WeightedMeasurement;
> import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;
> public class QuadraticFitterProblem {
>        private QuadraticFitter quadraticFitter;
>        public QuadraticFitterProblem() {
>          // create the uninitialized fitting problem
>          quadraticFitter = new QuadraticFitter();
>          quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
>          quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
>          quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
>          quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
>          quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
>          try {
>            // solve the problem, using a Levenberg-Marquardt algorithm with
> default settings
>            LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
>            //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
>            estimator.estimate(quadraticFitter);
>          } catch (EstimationException ee) {
>                System.err.println(ee.getMessage());
>          }
>        }
>        /**
>         * @param args
>         *
>         */
>        public static void main(String[] args) {
>                        new QuadraticFitterProblem();
>                        System.out.println ("Done.");
>        }
> }
> ----------------------------------------------------------------------------------------------
> import org.apache.commons.math.estimation.EstimatedParameter;
> //import org.apache.commons.math.estimation.EstimationException;
> //import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> import org.apache.commons.math.estimation.SimpleEstimationProblem;
> import org.apache.commons.math.estimation.WeightedMeasurement;
> public class QuadraticFitter extends SimpleEstimationProblem {
>        // y = a x<sup>2</sup> + b x + c
>    private EstimatedParameter a;
>    private EstimatedParameter b;
>    private EstimatedParameter c;
>    /**
>     * constructor
>     *
>     *Fitter for a quadratic model to a sample of 2D points.
>     * <p>The model is y(x) = a x<sup>2</sup> + b x + c
>     * its three parameters of the model are a, b and c.</p>
>     */
>    public QuadraticFitter() {
>        // three parameters of the model
>        a = new EstimatedParameter("a", 0.0);
>        b = new EstimatedParameter("b", 0.0);
>        c = new EstimatedParameter("c", 0.0);
>        // provide the parameters to the base class which
>        // implements the getAllParameters and getUnboundParameters methods
>        addParameter(a);
>        addParameter(b);
>        addParameter(c);
>    }
>    /**
>     * Add a sample point
>     *
>     * @param x abscissa
>     * @param y ordinate
>     * @param w weight
>     */
>    public void addPoint(double x, double y, double w) {
>        addMeasurement(new LocalMeasurement(x, y, w));
>    }
>    /**
>     * Get the value of the quadratic coefficient.
>     *
>     * @return the value of a for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getA() {
>        return a.getEstimate();
>    }
>    /**
>     * Get the value of the linear coefficient.
>     *
>     * @return the value of b for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getB() {
>        return b.getEstimate();
>    }
>    /**
>     * Get the value of the constant coefficient.
>     *
>     * @return the value of ac for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getC() {
>        return c.getEstimate();
>    }
>    /**
>     * Get the theoretical value of the model for some x.
>     * <p>The theoretical value is the value computed using
>     * the current state of the problem parameters.</p>
>     *
>     * Note the use of Hörner's method (synthetic division) for
> evaluating polynomials,
>     * (more efficient)
>     *
>     * @param x explanatory variable
>     * @return the theoretical value y = a x<sup>2</sup> + b x + c
>     */
>    public double theoreticalValue(double x) {
>        //System.out.println ("x = " + x + "  a.getEstimate() = " +
> a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
> c.getEstimate() = " + c.getEstimate());
>        return ( (a.getEstimate() * x + b.getEstimate() ) * x +
> c.getEstimate());
>    }
>    /**
>     * Get the partial derivative of the theoretical value
>     * of the model for some x.
>     * <p>The derivative is computed using
>     * the current state of the problem parameters.</p>
>     *
>     * @param x explanatory variable
>     * @param parameter estimated parameter (either a, b, or c)
>     * @return the partial derivative dy/dp
>     */
>    private double partial(double x, EstimatedParameter parameter) {
>        // since we know the only parameters are a, b and c in this
>        // class we simply use "==" for efficiency
>        if (parameter == a) {
>            return x * x;
>        } else if (parameter == b) {
>            return x;
>        } else {
>            return 1.0;
>        }
>    }
>    /** Internal measurements class.
>     * <p>The measurement is the y value for a fixed specified x.</p>
>     */
>    private class LocalMeasurement extends WeightedMeasurement {
>        static final long serialVersionUID = 1;
>        private final double x;
>        // constructor
>        public LocalMeasurement(double x, double y, double w) {
>            super(w, y);
>            this.x = x;
>        }
>        public double getTheoreticalValue() {
>            // the value is provided by the model for the local x
>            return theoreticalValue(x);
>        }
>        public double getPartial(EstimatedParameter parameter) {
>            // the value is provided by the model for the local x
>            return partial(x, parameter);
>        }
>    }
>  }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.


[jira] Commented: (MATH-199) exception in LevenbergMarquardtEstimator

Posted by "Luc Maisonobe (JIRA)" <ji...@apache.org>.
    [ https://issues.apache.org/jira/browse/MATH-199?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12581590#action_12581590 ] 

Luc Maisonobe commented on MATH-199:
------------------------------------

Numbers that behave well are "normal numbers" as defined by IEEE754. Numbers to avoid in your case are "subnormal numbers". You can use Double.MIN_NORMAL if you use Java 6, or simply 4 / Double.MAX_VALUE (which is almost but not exactly the same value) if you use an earlier version of Java.

Beware that this may not be enough because you can test only the numbers you put in the method, not the intermediate values computed within it. For example the ak2, akk and alpha values are not available in the statement above, and betak is available only on output.

I'm not sure such constants have a place in [math], but this is only a personal opinion.

> exception in LevenbergMarquardtEstimator
> ----------------------------------------
>
>                 Key: MATH-199
>                 URL: https://issues.apache.org/jira/browse/MATH-199
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 1.2
>         Environment: Windows XP
> Java 6
>            Reporter: Mick
>
> I get this exception:
> Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
>        at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
>        at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
> on the code below.
> The exception does not occur all the weights in the quadraticFitter are 0.0;
> ---------------------------------------------------------------------------------------------
> package quadraticFitterProblem;
> import org.apache.commons.math.estimation.EstimationException;
> import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> //import org.apache.commons.math.estimation.WeightedMeasurement;
> import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;
> public class QuadraticFitterProblem {
>        private QuadraticFitter quadraticFitter;
>        public QuadraticFitterProblem() {
>          // create the uninitialized fitting problem
>          quadraticFitter = new QuadraticFitter();
>          quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
>          quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
>          quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
>          quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
>          quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
>          try {
>            // solve the problem, using a Levenberg-Marquardt algorithm with
> default settings
>            LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
>            //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
>            estimator.estimate(quadraticFitter);
>          } catch (EstimationException ee) {
>                System.err.println(ee.getMessage());
>          }
>        }
>        /**
>         * @param args
>         *
>         */
>        public static void main(String[] args) {
>                        new QuadraticFitterProblem();
>                        System.out.println ("Done.");
>        }
> }
> ----------------------------------------------------------------------------------------------
> import org.apache.commons.math.estimation.EstimatedParameter;
> //import org.apache.commons.math.estimation.EstimationException;
> //import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> import org.apache.commons.math.estimation.SimpleEstimationProblem;
> import org.apache.commons.math.estimation.WeightedMeasurement;
> public class QuadraticFitter extends SimpleEstimationProblem {
>        // y = a x<sup>2</sup> + b x + c
>    private EstimatedParameter a;
>    private EstimatedParameter b;
>    private EstimatedParameter c;
>    /**
>     * constructor
>     *
>     *Fitter for a quadratic model to a sample of 2D points.
>     * <p>The model is y(x) = a x<sup>2</sup> + b x + c
>     * its three parameters of the model are a, b and c.</p>
>     */
>    public QuadraticFitter() {
>        // three parameters of the model
>        a = new EstimatedParameter("a", 0.0);
>        b = new EstimatedParameter("b", 0.0);
>        c = new EstimatedParameter("c", 0.0);
>        // provide the parameters to the base class which
>        // implements the getAllParameters and getUnboundParameters methods
>        addParameter(a);
>        addParameter(b);
>        addParameter(c);
>    }
>    /**
>     * Add a sample point
>     *
>     * @param x abscissa
>     * @param y ordinate
>     * @param w weight
>     */
>    public void addPoint(double x, double y, double w) {
>        addMeasurement(new LocalMeasurement(x, y, w));
>    }
>    /**
>     * Get the value of the quadratic coefficient.
>     *
>     * @return the value of a for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getA() {
>        return a.getEstimate();
>    }
>    /**
>     * Get the value of the linear coefficient.
>     *
>     * @return the value of b for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getB() {
>        return b.getEstimate();
>    }
>    /**
>     * Get the value of the constant coefficient.
>     *
>     * @return the value of ac for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getC() {
>        return c.getEstimate();
>    }
>    /**
>     * Get the theoretical value of the model for some x.
>     * <p>The theoretical value is the value computed using
>     * the current state of the problem parameters.</p>
>     *
>     * Note the use of Hörner's method (synthetic division) for
> evaluating polynomials,
>     * (more efficient)
>     *
>     * @param x explanatory variable
>     * @return the theoretical value y = a x<sup>2</sup> + b x + c
>     */
>    public double theoreticalValue(double x) {
>        //System.out.println ("x = " + x + "  a.getEstimate() = " +
> a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
> c.getEstimate() = " + c.getEstimate());
>        return ( (a.getEstimate() * x + b.getEstimate() ) * x +
> c.getEstimate());
>    }
>    /**
>     * Get the partial derivative of the theoretical value
>     * of the model for some x.
>     * <p>The derivative is computed using
>     * the current state of the problem parameters.</p>
>     *
>     * @param x explanatory variable
>     * @param parameter estimated parameter (either a, b, or c)
>     * @return the partial derivative dy/dp
>     */
>    private double partial(double x, EstimatedParameter parameter) {
>        // since we know the only parameters are a, b and c in this
>        // class we simply use "==" for efficiency
>        if (parameter == a) {
>            return x * x;
>        } else if (parameter == b) {
>            return x;
>        } else {
>            return 1.0;
>        }
>    }
>    /** Internal measurements class.
>     * <p>The measurement is the y value for a fixed specified x.</p>
>     */
>    private class LocalMeasurement extends WeightedMeasurement {
>        static final long serialVersionUID = 1;
>        private final double x;
>        // constructor
>        public LocalMeasurement(double x, double y, double w) {
>            super(w, y);
>            this.x = x;
>        }
>        public double getTheoreticalValue() {
>            // the value is provided by the model for the local x
>            return theoreticalValue(x);
>        }
>        public double getPartial(EstimatedParameter parameter) {
>            // the value is provided by the model for the local x
>            return partial(x, parameter);
>        }
>    }
>  }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.


[jira] Resolved: (MATH-199) exception in LevenbergMarquardtEstimator

Posted by "Luc Maisonobe (JIRA)" <ji...@apache.org>.
     [ https://issues.apache.org/jira/browse/MATH-199?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Luc Maisonobe resolved MATH-199.
--------------------------------

    Resolution: Fixed

Fixed in svn as of r640205

The problem was due to an overflow in Q.R decomposition. One of  the transformed columns had both infinite and NaN elements, so the test of the norm was never met and a column index was never set.

The fix consist in detecting non-numeric norms and throwing an EstimationException stating the Q.R decomposition could not be performed.

> exception in LevenbergMarquardtEstimator
> ----------------------------------------
>
>                 Key: MATH-199
>                 URL: https://issues.apache.org/jira/browse/MATH-199
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 1.2
>         Environment: Windows XP
> Java 6
>            Reporter: Mick
>
> I get this exception:
> Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
>        at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
>        at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
> on the code below.
> The exception does not occur all the weights in the quadraticFitter are 0.0;
> ---------------------------------------------------------------------------------------------
> package quadraticFitterProblem;
> import org.apache.commons.math.estimation.EstimationException;
> import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> //import org.apache.commons.math.estimation.WeightedMeasurement;
> import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;
> public class QuadraticFitterProblem {
>        private QuadraticFitter quadraticFitter;
>        public QuadraticFitterProblem() {
>          // create the uninitialized fitting problem
>          quadraticFitter = new QuadraticFitter();
>          quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
>          quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
>          quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
>          quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
>          quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
>          try {
>            // solve the problem, using a Levenberg-Marquardt algorithm with
> default settings
>            LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
>            //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
>            estimator.estimate(quadraticFitter);
>          } catch (EstimationException ee) {
>                System.err.println(ee.getMessage());
>          }
>        }
>        /**
>         * @param args
>         *
>         */
>        public static void main(String[] args) {
>                        new QuadraticFitterProblem();
>                        System.out.println ("Done.");
>        }
> }
> ----------------------------------------------------------------------------------------------
> import org.apache.commons.math.estimation.EstimatedParameter;
> //import org.apache.commons.math.estimation.EstimationException;
> //import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> import org.apache.commons.math.estimation.SimpleEstimationProblem;
> import org.apache.commons.math.estimation.WeightedMeasurement;
> public class QuadraticFitter extends SimpleEstimationProblem {
>        // y = a x<sup>2</sup> + b x + c
>    private EstimatedParameter a;
>    private EstimatedParameter b;
>    private EstimatedParameter c;
>    /**
>     * constructor
>     *
>     *Fitter for a quadratic model to a sample of 2D points.
>     * <p>The model is y(x) = a x<sup>2</sup> + b x + c
>     * its three parameters of the model are a, b and c.</p>
>     */
>    public QuadraticFitter() {
>        // three parameters of the model
>        a = new EstimatedParameter("a", 0.0);
>        b = new EstimatedParameter("b", 0.0);
>        c = new EstimatedParameter("c", 0.0);
>        // provide the parameters to the base class which
>        // implements the getAllParameters and getUnboundParameters methods
>        addParameter(a);
>        addParameter(b);
>        addParameter(c);
>    }
>    /**
>     * Add a sample point
>     *
>     * @param x abscissa
>     * @param y ordinate
>     * @param w weight
>     */
>    public void addPoint(double x, double y, double w) {
>        addMeasurement(new LocalMeasurement(x, y, w));
>    }
>    /**
>     * Get the value of the quadratic coefficient.
>     *
>     * @return the value of a for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getA() {
>        return a.getEstimate();
>    }
>    /**
>     * Get the value of the linear coefficient.
>     *
>     * @return the value of b for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getB() {
>        return b.getEstimate();
>    }
>    /**
>     * Get the value of the constant coefficient.
>     *
>     * @return the value of ac for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getC() {
>        return c.getEstimate();
>    }
>    /**
>     * Get the theoretical value of the model for some x.
>     * <p>The theoretical value is the value computed using
>     * the current state of the problem parameters.</p>
>     *
>     * Note the use of Hörner's method (synthetic division) for
> evaluating polynomials,
>     * (more efficient)
>     *
>     * @param x explanatory variable
>     * @return the theoretical value y = a x<sup>2</sup> + b x + c
>     */
>    public double theoreticalValue(double x) {
>        //System.out.println ("x = " + x + "  a.getEstimate() = " +
> a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
> c.getEstimate() = " + c.getEstimate());
>        return ( (a.getEstimate() * x + b.getEstimate() ) * x +
> c.getEstimate());
>    }
>    /**
>     * Get the partial derivative of the theoretical value
>     * of the model for some x.
>     * <p>The derivative is computed using
>     * the current state of the problem parameters.</p>
>     *
>     * @param x explanatory variable
>     * @param parameter estimated parameter (either a, b, or c)
>     * @return the partial derivative dy/dp
>     */
>    private double partial(double x, EstimatedParameter parameter) {
>        // since we know the only parameters are a, b and c in this
>        // class we simply use "==" for efficiency
>        if (parameter == a) {
>            return x * x;
>        } else if (parameter == b) {
>            return x;
>        } else {
>            return 1.0;
>        }
>    }
>    /** Internal measurements class.
>     * <p>The measurement is the y value for a fixed specified x.</p>
>     */
>    private class LocalMeasurement extends WeightedMeasurement {
>        static final long serialVersionUID = 1;
>        private final double x;
>        // constructor
>        public LocalMeasurement(double x, double y, double w) {
>            super(w, y);
>            this.x = x;
>        }
>        public double getTheoreticalValue() {
>            // the value is provided by the model for the local x
>            return theoreticalValue(x);
>        }
>        public double getPartial(EstimatedParameter parameter) {
>            // the value is provided by the model for the local x
>            return partial(x, parameter);
>        }
>    }
>  }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.


[jira] Commented: (MATH-199) exception in LevenbergMarquardtEstimator

Posted by "Luc Maisonobe (JIRA)" <ji...@apache.org>.
    [ https://issues.apache.org/jira/browse/MATH-199?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12581558#action_12581558 ] 

Luc Maisonobe commented on MATH-199:
------------------------------------

Yes, this is this small value that triggers the error. It is a valid number but during the Q.R decomposition, this number is used in several operations. One of these operations lead to an overflow.

At line 785 of LevenbergMarquardt.java, we compute :  double betak = 1.0 / (ak2 - akk * alpha)

With the very small values of the example, the value for betak overflows while processing the second column (first column is computed without problem). The denominator ak2 - akk * alpha has a value of 1.43e-322 which can be handled by double. It's inverse is 6.9e321, far larger than the Double.MAX_VALUE which is about 1.8e+308. The result of the computation is that betak is set to positive infinity. The rest of the computation behaves badly with such a value. Some of the elements of third column are set to infinity, others are set to NaN.

The fact is that despite many representable double number have a representable inverse in IEE754, it is not true for all.

> exception in LevenbergMarquardtEstimator
> ----------------------------------------
>
>                 Key: MATH-199
>                 URL: https://issues.apache.org/jira/browse/MATH-199
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 1.2
>         Environment: Windows XP
> Java 6
>            Reporter: Mick
>
> I get this exception:
> Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.qrDecomposition(LevenbergMarquardtEstimator.java:772)
>        at org.apache.commons.math.estimation.LevenbergMarquardtEstimator.estimate(LevenbergMarquardtEstimator.java:232)
>        at quadraticFitterProblem.QuadraticFitterProblem.<init>(QuadraticFitterProblem.java:27)
>        at quadraticFitterProblem.QuadraticFitterProblem.main(QuadraticFitterProblem.java:40)
> on the code below.
> The exception does not occur all the weights in the quadraticFitter are 0.0;
> ---------------------------------------------------------------------------------------------
> package quadraticFitterProblem;
> import org.apache.commons.math.estimation.EstimationException;
> import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> //import org.apache.commons.math.estimation.WeightedMeasurement;
> import com.strategicanalytics.dtd.data.smoothers.QuadraticFitter;
> public class QuadraticFitterProblem {
>        private QuadraticFitter quadraticFitter;
>        public QuadraticFitterProblem() {
>          // create the uninitialized fitting problem
>          quadraticFitter = new QuadraticFitter();
>          quadraticFitter.addPoint (0,  -3.182591015485607, 0.0);
>          quadraticFitter.addPoint (1,  -2.5581184967730577, 4.4E-323);
>          quadraticFitter.addPoint (2,  -2.1488478161387325, 1.0);
>          quadraticFitter.addPoint (3,  -1.9122489313410047, 4.4E-323);
>          quadraticFitter.addPoint (4,  1.7785661310051026, 0.0);
>          try {
>            // solve the problem, using a Levenberg-Marquardt algorithm with
> default settings
>            LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
>            //WeightedMeasurement[] wm = quadraticFitter.getMeasurements();
>            estimator.estimate(quadraticFitter);
>          } catch (EstimationException ee) {
>                System.err.println(ee.getMessage());
>          }
>        }
>        /**
>         * @param args
>         *
>         */
>        public static void main(String[] args) {
>                        new QuadraticFitterProblem();
>                        System.out.println ("Done.");
>        }
> }
> ----------------------------------------------------------------------------------------------
> import org.apache.commons.math.estimation.EstimatedParameter;
> //import org.apache.commons.math.estimation.EstimationException;
> //import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
> import org.apache.commons.math.estimation.SimpleEstimationProblem;
> import org.apache.commons.math.estimation.WeightedMeasurement;
> public class QuadraticFitter extends SimpleEstimationProblem {
>        // y = a x<sup>2</sup> + b x + c
>    private EstimatedParameter a;
>    private EstimatedParameter b;
>    private EstimatedParameter c;
>    /**
>     * constructor
>     *
>     *Fitter for a quadratic model to a sample of 2D points.
>     * <p>The model is y(x) = a x<sup>2</sup> + b x + c
>     * its three parameters of the model are a, b and c.</p>
>     */
>    public QuadraticFitter() {
>        // three parameters of the model
>        a = new EstimatedParameter("a", 0.0);
>        b = new EstimatedParameter("b", 0.0);
>        c = new EstimatedParameter("c", 0.0);
>        // provide the parameters to the base class which
>        // implements the getAllParameters and getUnboundParameters methods
>        addParameter(a);
>        addParameter(b);
>        addParameter(c);
>    }
>    /**
>     * Add a sample point
>     *
>     * @param x abscissa
>     * @param y ordinate
>     * @param w weight
>     */
>    public void addPoint(double x, double y, double w) {
>        addMeasurement(new LocalMeasurement(x, y, w));
>    }
>    /**
>     * Get the value of the quadratic coefficient.
>     *
>     * @return the value of a for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getA() {
>        return a.getEstimate();
>    }
>    /**
>     * Get the value of the linear coefficient.
>     *
>     * @return the value of b for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getB() {
>        return b.getEstimate();
>    }
>    /**
>     * Get the value of the constant coefficient.
>     *
>     * @return the value of ac for the quadratic model
>     * y = a x<sup>2</sup> + b x + c
>     */
>    public double getC() {
>        return c.getEstimate();
>    }
>    /**
>     * Get the theoretical value of the model for some x.
>     * <p>The theoretical value is the value computed using
>     * the current state of the problem parameters.</p>
>     *
>     * Note the use of Hörner's method (synthetic division) for
> evaluating polynomials,
>     * (more efficient)
>     *
>     * @param x explanatory variable
>     * @return the theoretical value y = a x<sup>2</sup> + b x + c
>     */
>    public double theoreticalValue(double x) {
>        //System.out.println ("x = " + x + "  a.getEstimate() = " +
> a.getEstimate() + "  b.getEstimate() = " + b.getEstimate() + "
> c.getEstimate() = " + c.getEstimate());
>        return ( (a.getEstimate() * x + b.getEstimate() ) * x +
> c.getEstimate());
>    }
>    /**
>     * Get the partial derivative of the theoretical value
>     * of the model for some x.
>     * <p>The derivative is computed using
>     * the current state of the problem parameters.</p>
>     *
>     * @param x explanatory variable
>     * @param parameter estimated parameter (either a, b, or c)
>     * @return the partial derivative dy/dp
>     */
>    private double partial(double x, EstimatedParameter parameter) {
>        // since we know the only parameters are a, b and c in this
>        // class we simply use "==" for efficiency
>        if (parameter == a) {
>            return x * x;
>        } else if (parameter == b) {
>            return x;
>        } else {
>            return 1.0;
>        }
>    }
>    /** Internal measurements class.
>     * <p>The measurement is the y value for a fixed specified x.</p>
>     */
>    private class LocalMeasurement extends WeightedMeasurement {
>        static final long serialVersionUID = 1;
>        private final double x;
>        // constructor
>        public LocalMeasurement(double x, double y, double w) {
>            super(w, y);
>            this.x = x;
>        }
>        public double getTheoreticalValue() {
>            // the value is provided by the model for the local x
>            return theoreticalValue(x);
>        }
>        public double getPartial(EstimatedParameter parameter) {
>            // the value is provided by the model for the local x
>            return partial(x, parameter);
>        }
>    }
>  }

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.