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   * @since 2.0
49   * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
50   * {@link WeightedObservedPoints} instead.
51   */
52  @Deprecated
53  public class CurveFitter<T extends ParametricUnivariateFunction> {
54      /** Optimizer to use for the fitting. */
55      private final MultivariateVectorOptimizer optimizer;
56      /** Observed points. */
57      private final List<WeightedObservedPoint> observations;
58  
59      /**
60       * Simple constructor.
61       *
62       * @param optimizer Optimizer to use for the fitting.
63       * @since 3.1
64       */
65      public CurveFitter(final MultivariateVectorOptimizer optimizer) {
66          this.optimizer = optimizer;
67          observations = new ArrayList<WeightedObservedPoint>();
68      }
69  
70      /** Add an observed (x,y) point to the sample with unit weight.
71       * <p>Calling this method is equivalent to call
72       * {@code addObservedPoint(1.0, x, y)}.</p>
73       * @param x abscissa of the point
74       * @param y observed value of the point at x, after fitting we should
75       * have f(x) as close as possible to this value
76       * @see #addObservedPoint(double, double, double)
77       * @see #addObservedPoint(WeightedObservedPoint)
78       * @see #getObservations()
79       */
80      public void addObservedPoint(double x, double y) {
81          addObservedPoint(1.0, x, y);
82      }
83  
84      /** Add an observed weighted (x,y) point to the sample.
85       * @param weight weight of the observed point in the fit
86       * @param x abscissa of the point
87       * @param y observed value of the point at x, after fitting we should
88       * have f(x) as close as possible to this value
89       * @see #addObservedPoint(double, double)
90       * @see #addObservedPoint(WeightedObservedPoint)
91       * @see #getObservations()
92       */
93      public void addObservedPoint(double weight, double x, double y) {
94          observations.add(new WeightedObservedPoint(weight, x, y));
95      }
96  
97      /** Add an observed weighted (x,y) point to the sample.
98       * @param observed observed point to add
99       * @see #addObservedPoint(double, double)
100      * @see #addObservedPoint(double, double, double)
101      * @see #getObservations()
102      */
103     public void addObservedPoint(WeightedObservedPoint observed) {
104         observations.add(observed);
105     }
106 
107     /** Get the observed points.
108      * @return observed points
109      * @see #addObservedPoint(double, double)
110      * @see #addObservedPoint(double, double, double)
111      * @see #addObservedPoint(WeightedObservedPoint)
112      */
113     public WeightedObservedPoint[] getObservations() {
114         return observations.toArray(new WeightedObservedPoint[observations.size()]);
115     }
116 
117     /**
118      * Remove all observations.
119      */
120     public void clearObservations() {
121         observations.clear();
122     }
123 
124     /**
125      * Fit a curve.
126      * This method compute the coefficients of the curve that best
127      * fit the sample of observed points previously given through calls
128      * to the {@link #addObservedPoint(WeightedObservedPoint)
129      * addObservedPoint} method.
130      *
131      * @param f parametric function to fit.
132      * @param initialGuess first guess of the function parameters.
133      * @return the fitted parameters.
134      * @throws org.apache.commons.math3.exception.DimensionMismatchException
135      * if the start point dimension is wrong.
136      */
137     public double[] fit(T f, final double[] initialGuess) {
138         return fit(Integer.MAX_VALUE, f, initialGuess);
139     }
140 
141     /**
142      * Fit a curve.
143      * This method compute the coefficients of the curve that best
144      * fit the sample of observed points previously given through calls
145      * to the {@link #addObservedPoint(WeightedObservedPoint)
146      * addObservedPoint} method.
147      *
148      * @param f parametric function to fit.
149      * @param initialGuess first guess of the function parameters.
150      * @param maxEval Maximum number of function evaluations.
151      * @return the fitted parameters.
152      * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
153      * if the number of allowed evaluations is exceeded.
154      * @throws org.apache.commons.math3.exception.DimensionMismatchException
155      * if the start point dimension is wrong.
156      * @since 3.0
157      */
158     public double[] fit(int maxEval, T f,
159                         final double[] initialGuess) {
160         // Prepare least squares problem.
161         double[] target  = new double[observations.size()];
162         double[] weights = new double[observations.size()];
163         int i = 0;
164         for (WeightedObservedPoint point : observations) {
165             target[i]  = point.getY();
166             weights[i] = point.getWeight();
167             ++i;
168         }
169 
170         // Input to the optimizer: the model and its Jacobian.
171         final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);
172 
173         // Perform the fit.
174         final PointVectorValuePair optimum
175             = optimizer.optimize(new MaxEval(maxEval),
176                                  model.getModelFunction(),
177                                  model.getModelFunctionJacobian(),
178                                  new Target(target),
179                                  new Weight(weights),
180                                  new InitialGuess(initialGuess));
181         // Extract the coefficients.
182         return optimum.getPointRef();
183     }
184 
185     /** Vectorial function computing function theoretical values. */
186     private class TheoreticalValuesFunction {
187         /** Function to fit. */
188         private final ParametricUnivariateFunction f;
189 
190         /**
191          * @param f function to fit.
192          */
193         public TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
194             this.f = f;
195         }
196 
197         /**
198          * @return the model function values.
199          */
200         public ModelFunction getModelFunction() {
201             return new ModelFunction(new MultivariateVectorFunction() {
202                     /** {@inheritDoc} */
203                     public double[] value(double[] point) {
204                         // compute the residuals
205                         final double[] values = new double[observations.size()];
206                         int i = 0;
207                         for (WeightedObservedPoint observed : observations) {
208                             values[i++] = f.value(observed.getX(), point);
209                         }
210 
211                         return values;
212                     }
213                 });
214         }
215 
216         /**
217          * @return the model function Jacobian.
218          */
219         public ModelFunctionJacobian getModelFunctionJacobian() {
220             return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
221                     public double[][] value(double[] point) {
222                         final double[][] jacobian = new double[observations.size()][];
223                         int i = 0;
224                         for (WeightedObservedPoint observed : observations) {
225                             jacobian[i++] = f.gradient(observed.getX(), point);
226                         }
227                         return jacobian;
228                     }
229                 });
230         }
231     }
232 }