View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.fitting;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  import org.apache.commons.math3.analysis.MultivariateVectorFunction;
22  import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
23  import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
24  import org.apache.commons.math3.optim.MaxEval;
25  import org.apache.commons.math3.optim.InitialGuess;
26  import org.apache.commons.math3.optim.PointVectorValuePair;
27  import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
28  import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
29  import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
30  import org.apache.commons.math3.optim.nonlinear.vector.Target;
31  import org.apache.commons.math3.optim.nonlinear.vector.Weight;
32  
33  /**
34   * Fitter for parametric univariate real functions y = f(x).
35   * <br/>
36   * When a univariate real function y = f(x) does depend on some
37   * unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
38   * this class can be used to find these parameters. It does this
39   * by <em>fitting</em> the curve so it remains very close to a set of
40   * observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
41   * y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
42   * is done by finding the parameters values that minimizes the objective
43   * function &sum;(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
44   * really a least squares problem.
45   *
46   * @param <T> Function to use for the fit.
47   *
48   * @version $Id: CurveFitter.java 1416643 2012-12-03 19:37:14Z tn $
49   * @since 2.0
50   * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
51   * {@link WeightedObservedPoints} instead.
52   */
53  @Deprecated
54  public class CurveFitter<T extends ParametricUnivariateFunction> {
55      /** Optimizer to use for the fitting. */
56      private final MultivariateVectorOptimizer optimizer;
57      /** Observed points. */
58      private final List<WeightedObservedPoint> observations;
59  
60      /**
61       * Simple constructor.
62       *
63       * @param optimizer Optimizer to use for the fitting.
64       * @since 3.1
65       */
66      public CurveFitter(final MultivariateVectorOptimizer optimizer) {
67          this.optimizer = optimizer;
68          observations = new ArrayList<WeightedObservedPoint>();
69      }
70  
71      /** Add an observed (x,y) point to the sample with unit weight.
72       * <p>Calling this method is equivalent to call
73       * {@code addObservedPoint(1.0, x, y)}.</p>
74       * @param x abscissa of the point
75       * @param y observed value of the point at x, after fitting we should
76       * have f(x) as close as possible to this value
77       * @see #addObservedPoint(double, double, double)
78       * @see #addObservedPoint(WeightedObservedPoint)
79       * @see #getObservations()
80       */
81      public void addObservedPoint(double x, double y) {
82          addObservedPoint(1.0, x, y);
83      }
84  
85      /** Add an observed weighted (x,y) point to the sample.
86       * @param weight weight of the observed point in the fit
87       * @param x abscissa of the point
88       * @param y observed value of the point at x, after fitting we should
89       * have f(x) as close as possible to this value
90       * @see #addObservedPoint(double, double)
91       * @see #addObservedPoint(WeightedObservedPoint)
92       * @see #getObservations()
93       */
94      public void addObservedPoint(double weight, double x, double y) {
95          observations.add(new WeightedObservedPoint(weight, x, y));
96      }
97  
98      /** Add an observed weighted (x,y) point to the sample.
99       * @param observed observed point to add
100      * @see #addObservedPoint(double, double)
101      * @see #addObservedPoint(double, double, double)
102      * @see #getObservations()
103      */
104     public void addObservedPoint(WeightedObservedPoint observed) {
105         observations.add(observed);
106     }
107 
108     /** Get the observed points.
109      * @return observed points
110      * @see #addObservedPoint(double, double)
111      * @see #addObservedPoint(double, double, double)
112      * @see #addObservedPoint(WeightedObservedPoint)
113      */
114     public WeightedObservedPoint[] getObservations() {
115         return observations.toArray(new WeightedObservedPoint[observations.size()]);
116     }
117 
118     /**
119      * Remove all observations.
120      */
121     public void clearObservations() {
122         observations.clear();
123     }
124 
125     /**
126      * Fit a curve.
127      * This method compute the coefficients of the curve that best
128      * fit the sample of observed points previously given through calls
129      * to the {@link #addObservedPoint(WeightedObservedPoint)
130      * addObservedPoint} method.
131      *
132      * @param f parametric function to fit.
133      * @param initialGuess first guess of the function parameters.
134      * @return the fitted parameters.
135      * @throws org.apache.commons.math3.exception.DimensionMismatchException
136      * if the start point dimension is wrong.
137      */
138     public double[] fit(T f, final double[] initialGuess) {
139         return fit(Integer.MAX_VALUE, f, initialGuess);
140     }
141 
142     /**
143      * Fit a curve.
144      * This method compute the coefficients of the curve that best
145      * fit the sample of observed points previously given through calls
146      * to the {@link #addObservedPoint(WeightedObservedPoint)
147      * addObservedPoint} method.
148      *
149      * @param f parametric function to fit.
150      * @param initialGuess first guess of the function parameters.
151      * @param maxEval Maximum number of function evaluations.
152      * @return the fitted parameters.
153      * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
154      * if the number of allowed evaluations is exceeded.
155      * @throws org.apache.commons.math3.exception.DimensionMismatchException
156      * if the start point dimension is wrong.
157      * @since 3.0
158      */
159     public double[] fit(int maxEval, T f,
160                         final double[] initialGuess) {
161         // Prepare least squares problem.
162         double[] target  = new double[observations.size()];
163         double[] weights = new double[observations.size()];
164         int i = 0;
165         for (WeightedObservedPoint point : observations) {
166             target[i]  = point.getY();
167             weights[i] = point.getWeight();
168             ++i;
169         }
170 
171         // Input to the optimizer: the model and its Jacobian.
172         final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);
173 
174         // Perform the fit.
175         final PointVectorValuePair optimum
176             = optimizer.optimize(new MaxEval(maxEval),
177                                  model.getModelFunction(),
178                                  model.getModelFunctionJacobian(),
179                                  new Target(target),
180                                  new Weight(weights),
181                                  new InitialGuess(initialGuess));
182         // Extract the coefficients.
183         return optimum.getPointRef();
184     }
185 
186     /** Vectorial function computing function theoretical values. */
187     private class TheoreticalValuesFunction {
188         /** Function to fit. */
189         private final ParametricUnivariateFunction f;
190 
191         /**
192          * @param f function to fit.
193          */
194         public TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
195             this.f = f;
196         }
197 
198         /**
199          * @return the model function values.
200          */
201         public ModelFunction getModelFunction() {
202             return new ModelFunction(new MultivariateVectorFunction() {
203                     /** {@inheritDoc} */
204                     public double[] value(double[] point) {
205                         // compute the residuals
206                         final double[] values = new double[observations.size()];
207                         int i = 0;
208                         for (WeightedObservedPoint observed : observations) {
209                             values[i++] = f.value(observed.getX(), point);
210                         }
211 
212                         return values;
213                     }
214                 });
215         }
216 
217         /**
218          * @return the model function Jacobian.
219          */
220         public ModelFunctionJacobian getModelFunctionJacobian() {
221             return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
222                     public double[][] value(double[] point) {
223                         final double[][] jacobian = new double[observations.size()][];
224                         int i = 0;
225                         for (WeightedObservedPoint observed : observations) {
226                             jacobian[i++] = f.gradient(observed.getX(), point);
227                         }
228                         return jacobian;
229                     }
230                 });
231         }
232     }
233 }