7

How do I use Math Commons CurveFitter to fit a function to a set of data? I was told to use CurveFitter with LevenbergMarquardtOptimizer and ParametricUnivariateFunction, but I don't know what to write in the ParametricUnivariateFunction gradient and value methods. Besides, after writing them, how to get the fitted function parameters? My function:

public static double fnc(double t, double a, double b, double c){
  return a * Math.pow(t, b) * Math.exp(-c * t);
}

2 Answers 2

17

So, this is an old question, but I ran into the same issue recently, and ended up having to delve into mailing lists and the Apache Commons Math source code to figure it out.

This API is remarkably poorly documented, but in the current version of Apache Common Math (3.3+), there are two parts, assuming you have a single variable with multiple parameters: the function to fit with (which implements ParametricUnivariateFunction) and the curve fitter (which extends AbstractCurveFitter).

Function to Fit

  • public double value(double t, double... parameters)
    • Your equation. This is where you would put your fnc logic.
  • public double[] gradient(double t, double... parameters)
    • Returns an array of partial derivative of the above with respect to each parameters. This calculator may be helpful if (like me) you're rusty on your calculus, but any good computer algebra system can calculate these values.

Curve Fitter

  • protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> points)
    • Sets up a bunch of boilerplate crap, and returns a least squares problem for the fitter to use.

Putting it all together, here's an example solution in your specific case:

import java.util.*;
import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
import org.apache.commons.math3.fitting.AbstractCurveFitter;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem;
import org.apache.commons.math3.fitting.WeightedObservedPoint;
import org.apache.commons.math3.linear.DiagonalMatrix;

class MyFunc implements ParametricUnivariateFunction {
    public double value(double t, double... parameters) {
        return parameters[0] * Math.pow(t, parameters[1]) * Math.exp(-parameters[2] * t);
    }

    // Jacobian matrix of the above. In this case, this is just an array of
    // partial derivatives of the above function, with one element for each parameter.
    public double[] gradient(double t, double... parameters) {
        final double a = parameters[0];
        final double b = parameters[1];
        final double c = parameters[2];

        return new double[] {
            Math.exp(-c*t) * Math.pow(t, b),
            a * Math.exp(-c*t) * Math.pow(t, b) * Math.log(t),
            a * (-Math.exp(-c*t)) * Math.pow(t, b+1)
        };
    }
}

public class MyFuncFitter extends AbstractCurveFitter {
    protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> points) {
        final int len = points.size();
        final double[] target  = new double[len];
        final double[] weights = new double[len];
        final double[] initialGuess = { 1.0, 1.0, 1.0 };

        int i = 0;
        for(WeightedObservedPoint point : points) {
            target[i]  = point.getY();
            weights[i] = point.getWeight();
            i += 1;
        }

        final AbstractCurveFitter.TheoreticalValuesFunction model = new
            AbstractCurveFitter.TheoreticalValuesFunction(new MyFunc(), points);

        return new LeastSquaresBuilder().
            maxEvaluations(Integer.MAX_VALUE).
            maxIterations(Integer.MAX_VALUE).
            start(initialGuess).
            target(target).
            weight(new DiagonalMatrix(weights)).
            model(model.getModelFunction(), model.getModelFunctionJacobian()).
            build();
    }

    public static void main(String[] args) {
        MyFuncFitter fitter = new MyFuncFitter();
        ArrayList<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>();

        // Add points here; for instance,
        WeightedObservedPoint point = new WeightedObservedPoint(1.0,
            1.0,
            1.0);
        points.add(point);

        final double coeffs[] = fitter.fit(points);
        System.out.println(Arrays.toString(coeffs));
    }
}
2
  • 1
    I'm confused about the points collection. Is there no X_value to them?? Why does the target only contain the Y-value? Oct 12, 2015 at 15:45
  • Also how would I add constraints on parameters (for example that parameter a in f(x) = cln(ax) needs to always be positive ) ? Oct 14, 2015 at 15:57
5

I know this question's pretty old and i80and did an excellent job answering this, but I just thought to add (for future SO-ers) that there's a pretty easy way to compute derivatives or partial derivatives with Apache Math (so you don't have to do your own differentiation for the Jacobian Matrix). It's the DerivativeStructure.

Extending i80and's answer to use the DerivativeStructure class:

//Everything stays the same except for the Jacobian Matrix

import java.util.*;
import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
import org.apache.commons.math3.fitting.AbstractCurveFitter;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem;
import org.apache.commons.math3.fitting.WeightedObservedPoint;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;

class MyFunc implements ParametricUnivariateFunction {
    public double value(double t, double... parameters) {
        return parameters[0] * Math.pow(t, parameters[1]) * Math.exp(-parameters[2] * t);
    }

    // Jacobian matrix of the above. In this case, this is just an array of
    // partial derivatives of the above function, with one element for each parameter.
    public double[] gradient(double t, double... parameters) {
        final double a = parameters[0];
        final double b = parameters[1];
        final double c = parameters[2];

        // Jacobian Matrix Edit

        // Using Derivative Structures...
        // constructor takes 4 arguments - the number of parameters in your
        // equation to be differentiated (3 in this case), the order of
        // differentiation for the DerivativeStructure, the index of the
        // parameter represented by the DS, and the value of the parameter itself
        DerivativeStructure aDev = new DerivativeStructure(3, 1, 0, a);
        DerivativeStructure bDev = new DerivativeStructure(3, 1, 1, b);
        DerivativeStructure cDev = new DerivativeStructure(3, 1, 2, c);

        // define the equation to be differentiated using another DerivativeStructure
        DerivativeStructure y = aDev.multiply(DerivativeStructure.pow(t, bDev))
                .multiply(cDev.negate().multiply(t).exp());

        // then return the partial derivatives required
        // notice the format, 3 arguments for the method since 3 parameters were
        // specified first order derivative of the first parameter, then the second, 
        // then the third
        return new double[] {
                y.getPartialDerivative(1, 0, 0),
                y.getPartialDerivative(0, 1, 0),
                y.getPartialDerivative(0, 0, 1)
        };

    }
}

public class MyFuncFitter extends AbstractCurveFitter {
    protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> points) {
        final int len = points.size();
        final double[] target  = new double[len];
        final double[] weights = new double[len];
        final double[] initialGuess = { 1.0, 1.0, 1.0 };

        int i = 0;
        for(WeightedObservedPoint point : points) {
            target[i]  = point.getY();
            weights[i] = point.getWeight();
            i += 1;
        }

        final AbstractCurveFitter.TheoreticalValuesFunction model = new
                AbstractCurveFitter.TheoreticalValuesFunction(new MyFunc(), points);

        return new LeastSquaresBuilder().
                maxEvaluations(Integer.MAX_VALUE).
                maxIterations(Integer.MAX_VALUE).
                start(initialGuess).
                target(target).
                weight(new DiagonalMatrix(weights)).
                model(model.getModelFunction(), model.getModelFunctionJacobian()).
                build();
    }

    public static void main(String[] args) {
        MyFuncFitter fitter = new MyFuncFitter();
        ArrayList<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>();

        // Add points here; for instance,
        WeightedObservedPoint point = new WeightedObservedPoint(1.0,
                1.0,
                1.0);
        points.add(point);

        final double coeffs[] = fitter.fit(points);
        System.out.println(Arrays.toString(coeffs));
    }
}

And that's it. I know it's a pretty convoluted/confusing class to use, but it definitely comes in handy when you're dealing with very complicated equations that would be troublesome to get partial derivatives of by hand (this happened to me not so long ago), or when you want to derive partial derivatives say to the second or third order.

In the case of second, third, et cetera order derivatives, all you'll have to do is:

// specify the required order as the second argument, say second order so 2
DerivativeStructure aDev = new DerivativeStructure(3, 2, 0, a);        
DerivativeStructure bDev = new DerivativeStructure(3, 2, 1, b);
DerivativeStructure cDev = new DerivativeStructure(3, 2, 2, c);

// and then specify the order again here
y.getPartialDerivative(2, 0, 0),
y.getPartialDerivative(0, 2, 0),
y.getPartialDerivative(0, 0, 2)

Hopefully, this helps somebody sometime.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.