org.apache.commons.math3.optimization.fitting.HarmonicFitter Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of commons-math3 Show documentation
Show all versions of commons-math3 Show documentation
The Apache Commons Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang.
The newest version!
/*
* 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.
*/
package org.apache.commons.math3.optimization.fitting;
import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
import org.apache.commons.math3.analysis.function.HarmonicOscillator;
import org.apache.commons.math3.exception.ZeroException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.MathIllegalStateException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
/**
* Class that implements a curve fitting specialized for sinusoids.
*
* Harmonic fitting is a very simple case of curve fitting. The
* estimated coefficients are the amplitude a, the pulsation ω and
* the phase φ: f (t) = a cos (ω t + φ). They are
* searched by a least square estimator initialized with a rough guess
* based on integrals.
*
* @deprecated As of 3.1 (to be removed in 4.0).
* @since 2.0
*/
@Deprecated
public class HarmonicFitter extends CurveFitter {
/**
* Simple constructor.
* @param optimizer Optimizer to use for the fitting.
*/
public HarmonicFitter(final DifferentiableMultivariateVectorOptimizer optimizer) {
super(optimizer);
}
/**
* Fit an harmonic function to the observed points.
*
* @param initialGuess First guess values in the following order:
*
* - Amplitude
* - Angular frequency
* - Phase
*
* @return the parameters of the harmonic function that best fits the
* observed points (in the same order as above).
*/
public double[] fit(double[] initialGuess) {
return fit(new HarmonicOscillator.Parametric(), initialGuess);
}
/**
* Fit an harmonic function to the observed points.
* An initial guess will be automatically computed.
*
* @return the parameters of the harmonic function that best fits the
* observed points (see the other {@link #fit(double[]) fit} method.
* @throws NumberIsTooSmallException if the sample is too short for the
* the first guess to be computed.
* @throws ZeroException if the first guess cannot be computed because
* the abscissa range is zero.
*/
public double[] fit() {
return fit((new ParameterGuesser(getObservations())).guess());
}
/**
* This class guesses harmonic coefficients from a sample.
* The algorithm used to guess the coefficients is as follows:
*
* We know f (t) at some sampling points ti and want to find a,
* ω and φ such that f (t) = a cos (ω t + φ).
*
*
* From the analytical expression, we can compute two primitives :
*
* If2 (t) = ∫ f2 = a2 × [t + S (t)] / 2
* If'2 (t) = ∫ f'2 = a2 ω2 × [t - S (t)] / 2
* where S (t) = sin (2 (ω t + φ)) / (2 ω)
*
*
*
* We can remove S between these expressions :
*
* If'2 (t) = a2 ω2 t - ω2 If2 (t)
*
*
*
* The preceding expression shows that If'2 (t) is a linear
* combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t)
*
*
* From the primitive, we can deduce the same form for definite
* integrals between t1 and ti for each ti :
*
* If2 (ti) - If2 (t1) = A × (ti - t1) + B × (If2 (ti) - If2 (t1))
*
*
*
* We can find the coefficients A and B that best fit the sample
* to this linear expression by computing the definite integrals for
* each sample points.
*
*
* For a bilinear expression z (xi, yi) = A × xi + B × yi, the
* coefficients A and B that minimize a least square criterion
* ∑ (zi - z (xi, yi))2 are given by these expressions:
*
*
* ∑yiyi ∑xizi - ∑xiyi ∑yizi
* A = ------------------------
* ∑xixi ∑yiyi - ∑xiyi ∑xiyi
*
* ∑xixi ∑yizi - ∑xiyi ∑xizi
* B = ------------------------
* ∑xixi ∑yiyi - ∑xiyi ∑xiyi
*
*
*
*
* In fact, we can assume both a and ω are positive and
* compute them directly, knowing that A = a2 ω2 and that
* B = - ω2. The complete algorithm is therefore:
*
*
* for each ti from t1 to tn-1, compute:
* f (ti)
* f' (ti) = (f (ti+1) - f(ti-1)) / (ti+1 - ti-1)
* xi = ti - t1
* yi = ∫ f2 from t1 to ti
* zi = ∫ f'2 from t1 to ti
* update the sums ∑xixi, ∑yiyi, ∑xiyi, ∑xizi and ∑yizi
* end for
*
* |--------------------------
* \ | ∑yiyi ∑xizi - ∑xiyi ∑yizi
* a = \ | ------------------------
* \| ∑xiyi ∑xizi - ∑xixi ∑yizi
*
*
* |--------------------------
* \ | ∑xiyi ∑xizi - ∑xixi ∑yizi
* ω = \ | ------------------------
* \| ∑xixi ∑yiyi - ∑xiyi ∑xiyi
*
*
*
*
* Once we know ω, we can compute:
*
* fc = ω f (t) cos (ω t) - f' (t) sin (ω t)
* fs = ω f (t) sin (ω t) + f' (t) cos (ω t)
*
*
*
* It appears that fc = a ω cos (φ) and
* fs = -a ω sin (φ), so we can use these
* expressions to compute φ. The best estimate over the sample is
* given by averaging these expressions.
*
*
* Since integrals and means are involved in the preceding
* estimations, these operations run in O(n) time, where n is the
* number of measurements.
*/
public static class ParameterGuesser {
/** Amplitude. */
private final double a;
/** Angular frequency. */
private final double omega;
/** Phase. */
private final double phi;
/**
* Simple constructor.
*
* @param observations Sampled observations.
* @throws NumberIsTooSmallException if the sample is too short.
* @throws ZeroException if the abscissa range is zero.
* @throws MathIllegalStateException when the guessing procedure cannot
* produce sensible results.
*/
public ParameterGuesser(WeightedObservedPoint[] observations) {
if (observations.length < 4) {
throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
observations.length, 4, true);
}
final WeightedObservedPoint[] sorted = sortObservations(observations);
final double aOmega[] = guessAOmega(sorted);
a = aOmega[0];
omega = aOmega[1];
phi = guessPhi(sorted);
}
/**
* Gets an estimation of the parameters.
*
* @return the guessed parameters, in the following order:
*
* - Amplitude
* - Angular frequency
* - Phase
*
*/
public double[] guess() {
return new double[] { a, omega, phi };
}
/**
* Sort the observations with respect to the abscissa.
*
* @param unsorted Input observations.
* @return the input observations, sorted.
*/
private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
final WeightedObservedPoint[] observations = unsorted.clone();
// Since the samples are almost always already sorted, this
// method is implemented as an insertion sort that reorders the
// elements in place. Insertion sort is very efficient in this case.
WeightedObservedPoint curr = observations[0];
for (int j = 1; j < observations.length; ++j) {
WeightedObservedPoint prec = curr;
curr = observations[j];
if (curr.getX() < prec.getX()) {
// the current element should be inserted closer to the beginning
int i = j - 1;
WeightedObservedPoint mI = observations[i];
while ((i >= 0) && (curr.getX() < mI.getX())) {
observations[i + 1] = mI;
if (i-- != 0) {
mI = observations[i];
}
}
observations[i + 1] = curr;
curr = observations[j];
}
}
return observations;
}
/**
* Estimate a first guess of the amplitude and angular frequency.
* This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method
* has been called previously.
*
* @param observations Observations, sorted w.r.t. abscissa.
* @throws ZeroException if the abscissa range is zero.
* @throws MathIllegalStateException when the guessing procedure cannot
* produce sensible results.
* @return the guessed amplitude (at index 0) and circular frequency
* (at index 1).
*/
private double[] guessAOmega(WeightedObservedPoint[] observations) {
final double[] aOmega = new double[2];
// initialize the sums for the linear model between the two integrals
double sx2 = 0;
double sy2 = 0;
double sxy = 0;
double sxz = 0;
double syz = 0;
double currentX = observations[0].getX();
double currentY = observations[0].getY();
double f2Integral = 0;
double fPrime2Integral = 0;
final double startX = currentX;
for (int i = 1; i < observations.length; ++i) {
// one step forward
final double previousX = currentX;
final double previousY = currentY;
currentX = observations[i].getX();
currentY = observations[i].getY();
// update the integrals of f2 and f'2
// considering a linear model for f (and therefore constant f')
final double dx = currentX - previousX;
final double dy = currentY - previousY;
final double f2StepIntegral =
dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
final double fPrime2StepIntegral = dy * dy / dx;
final double x = currentX - startX;
f2Integral += f2StepIntegral;
fPrime2Integral += fPrime2StepIntegral;
sx2 += x * x;
sy2 += f2Integral * f2Integral;
sxy += x * f2Integral;
sxz += x * fPrime2Integral;
syz += f2Integral * fPrime2Integral;
}
// compute the amplitude and pulsation coefficients
double c1 = sy2 * sxz - sxy * syz;
double c2 = sxy * sxz - sx2 * syz;
double c3 = sx2 * sy2 - sxy * sxy;
if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
final int last = observations.length - 1;
// Range of the observations, assuming that the
// observations are sorted.
final double xRange = observations[last].getX() - observations[0].getX();
if (xRange == 0) {
throw new ZeroException();
}
aOmega[1] = 2 * Math.PI / xRange;
double yMin = Double.POSITIVE_INFINITY;
double yMax = Double.NEGATIVE_INFINITY;
for (int i = 1; i < observations.length; ++i) {
final double y = observations[i].getY();
if (y < yMin) {
yMin = y;
}
if (y > yMax) {
yMax = y;
}
}
aOmega[0] = 0.5 * (yMax - yMin);
} else {
if (c2 == 0) {
// In some ill-conditioned cases (cf. MATH-844), the guesser
// procedure cannot produce sensible results.
throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
}
aOmega[0] = FastMath.sqrt(c1 / c2);
aOmega[1] = FastMath.sqrt(c2 / c3);
}
return aOmega;
}
/**
* Estimate a first guess of the phase.
*
* @param observations Observations, sorted w.r.t. abscissa.
* @return the guessed phase.
*/
private double guessPhi(WeightedObservedPoint[] observations) {
// initialize the means
double fcMean = 0;
double fsMean = 0;
double currentX = observations[0].getX();
double currentY = observations[0].getY();
for (int i = 1; i < observations.length; ++i) {
// one step forward
final double previousX = currentX;
final double previousY = currentY;
currentX = observations[i].getX();
currentY = observations[i].getY();
final double currentYPrime = (currentY - previousY) / (currentX - previousX);
double omegaX = omega * currentX;
double cosine = FastMath.cos(omegaX);
double sine = FastMath.sin(omegaX);
fcMean += omega * currentY * cosine - currentYPrime * sine;
fsMean += omega * currentY * sine + currentYPrime * cosine;
}
return FastMath.atan2(-fsMean, fcMean);
}
}
}