Common non-linear curve fitting methods + more params

[Daniel15]

I used Fit.Curve for a four-parameter logistic (4PL) curve and it worked great

const double a = 0.0;
double[] xValues = new[] { 0.0, 57.0, 69.0, 80.0, 88.0, 97.0 };
double[] yValues = new[] { 0.0, 77.12, 88.78, 94.53, 96.83, 98.3 };

const double bGuess = 5;
const double cGuess = 43;
const double dGuess = 100;

var (b, c, d) = Fit.Curve(
    xValues,
    yValues,
    (b, c, d, x) => FourParameterLogisticCurve(a, b, c, d, x),
    initialGuess0: bGuess,
    initialGuess1: cGuess,
    initialGuess2: dGuess,
    maxIterations: 1000000
);
Console.WriteLine($"b={b}, c={c}, d={d}");

...
// https://www.myassays.com/four-parameter-logistic-regression.html
static double FourParameterLogisticCurve(double a, double b, double c, double d, double x)
{
    return d + ((a - d) / (1 + Math.Pow(x / c, b)));
}

Just some points of feedback:

• Would you be OK adding four parameter and five parameter overloads for Curve.Fit? Currently it only goes to three parameters, so I had to hard-code one of the parameter values. I can submit a PR for this if you like.
• What do you think about adding common curve functions to the library, so for example Fit.FourParameterLogistic, Fit.FiveParameterLogistic, etc

Thanks!

[Daniel15]

This is what I've got in one of my projects - Feel free to take any parts of it:

using System;
using MathNet.Numerics;
using System.Math;

namespace AtarCalc.Web.Utils
{
    /// <summary>
    /// Curve fitting algorithms
    /// </summary>
    public static class CurveFit
    {
        /// <summary>
        /// Four-Parameter Logistic (4PL) curve.
        /// Reference: https://www.myassays.com/four-parameter-logistic-regression.html
        /// </summary>
        public static double FourParameterLogisticCurve(double a, double b, double c, double d, double x)
        {
            return d + ((a - d) / (1 + System.Math.Pow(x / c, b)));
        }

        /// <summary>
        /// Four-Parameter Logistic (4PL) curve.
        /// Reference: https://www.myassays.com/four-parameter-logistic-regression.html
        /// </summary>
        public static double FourParameterLogisticCurve(FourParameterLogisticParams curveParams, double x)
        {
            return FourParameterLogisticCurve(curveParams.A, curveParams.B, curveParams.C, curveParams.D, x);
        }

        /// <summary>
        /// Computes the parameters for a Four-Parameter Logistic (4PL) curve.
        /// </summary>
        public static FourParameterLogisticParams FitFourParameterLogisticCurve(
            double[] xValues,
            double[] yValues,
            double bGuess,
            double cGuess,
            double dGuess,
            double a = 0
        )
        {
            var (b, c, d) = Fit.Curve(
                xValues,
                yValues,
                (b, c, d, x) => FourParameterLogisticCurve(a, b, c, d, x),
                initialGuess0: bGuess,
                initialGuess1: cGuess,
                initialGuess2: dGuess,
                maxIterations: 1000000
            );
            return new FourParameterLogisticParams(a, b, c, d);
        }

        /// <summary>
        /// Bell curve.
        /// y = a*e^(-(x - b)^2/(2*c^2))
        /// Reference: https://en.wikipedia.org/wiki/Bell_shaped_function
        /// </summary>
        /// <returns></returns>
        public static double GaussianBellCurve(double a, double b, double c, double x)
        {
            //
            return a * System.Math.Pow(System.Math.E, -(System.Math.Pow(x - b, 2) / (2 * System.Math.Pow(c, 2))));
        }

        /// <summary>
        /// Bell curve.
        /// y = a*e^(-(x - b)^2/(2*c^2))
        /// Reference: https://en.wikipedia.org/wiki/Bell_shaped_function
        /// </summary>
        /// <returns></returns>
        public static double GaussianBellCurve(BellCurveParams curveParams, double x)
        {
            return GaussianBellCurve(curveParams.A, curveParams.B, curveParams.C, x);
        }

        /// <summary>
        /// Computes the parameters for a Gaussian bell curve.
        /// </summary>
        public static BellCurveParams FitGaussianBellCurve(
            double[] xValues,
            double[] yValues,
            double aGuess,
            double bGuess,
            double cGuess
        )
        {
            var (a, b, c) = Fit.Curve(
                xValues,
                yValues,
                GaussianBellCurve,
                initialGuess0: aGuess,
                initialGuess1: bGuess,
                initialGuess2: cGuess,
                maxIterations: 1000000
            );
            return new BellCurveParams(a, b, c);
        }

        /// <summary>
        /// Power Curve
        /// y = ax^b
        /// </summary>
        public static double PowerCurve(double a, double b, double x)
        {
            return a * System.Math.Pow(x, b);
        }

        /// <summary>
        /// Power Curve
        /// y = ax^b
        /// </summary>
        public static double PowerCurve(PowerCurveParams curveParams, double x)
        {
            return PowerCurve(curveParams.A, curveParams.B, x);
        }

        /// <summary>
        /// Computes the parameters for a Power curve
        /// </summary>
        public static PowerCurveParams FitPowerCurve(
            double[] xValues,
            double[] yValues,
            double aGuess,
            double bGuess
        )
        {
            var (a, b) = Fit.Curve(
                xValues,
                yValues,
                PowerCurve,
                initialGuess0: aGuess,
                initialGuess1: bGuess,
                maxIterations: 1000000
            );
            return new PowerCurveParams(a, b);
        }

        /// <summary>
        /// Exponential - proportional growth
        /// y = Y0 - (V0/K)*(1 - e^(-K*x)
        /// </summary>
        public static double ExponentialProportionalGrowth(double y0, double v0, double k, double x)
        {
            return y0 - ((v0 / k) * (1 - System.Math.Pow(System.Math.E, -k * x)));
        }

        /// <summary>
        /// Exponential - proportional growth
        /// y = Y0 - (V0/K)*(1 - e^(-K*x)
        /// </summary>
        public static double ExponentialProportionalGrowth(ExponentialProportionalGrowthParams curveParams, double x)
        {
            return ExponentialProportionalGrowth(curveParams.Y0, curveParams.V0, curveParams.K, x);
        }

        /// <summary>
        /// Computes the parameters for an exponential proportional growth curve
        /// </summary>
        public static ExponentialProportionalGrowthParams FitExponentialProportionalGrowth(
            double[] xValues,
            double[] yValues,
            double y0Guess,
            double v0Guess,
            double kGuess
        )
        {
            var (y0, v0, k) = Fit.Curve(
                xValues,
                yValues,
                ExponentialProportionalGrowth,
                initialGuess0: y0Guess,
                initialGuess1: v0Guess,
                initialGuess2: kGuess,
                maxIterations: 1000000
            );
            return new ExponentialProportionalGrowthParams(y0, v0, k);
        }
    }

    public record FourParameterLogisticParams(double A, double B, double C, double D);

    public record BellCurveParams(double A, double B, double C);

    public record PowerCurveParams(double A, double B);

    public record ExponentialProportionalGrowthParams(double Y0, double V0, double K);
}

I return the parameters as a record so I can cache them easily. Given some particular source data, the app that uses this code computes the curve fit params lazily the first time they're needed in the app, then subsequent usages just get them from the cache rather than computing them again.

[Seebadori]

Hey @Daniel15 , Thx for sharing this.

I tried to use it with some of my assayData and i have error about 'Maximum iterations (1000000) reached.'

My values : double[] xValues = new[] { 5.85,13.8,28.5,66.0,124 };

double[] yValues = new[] { 2449.0,2131.0,1592.3,934.0,454 };

Can you give me some info why it can resolve this?

Actually i use R to resolve this. (with library(dr4pl)) but it will be perfect to use c# without R.Net

Thx and sorry for my english

[Daniel15]

'Maximum iterations (1000000) reached.'

@JulienHDev You may need better guess values. I'm not actually very familiar with this type of curve fitting method. Someone else here may have better knowledge than I do.

The other thing is that when I wrote that code, this library only supported three variables, so I had to hard-code one of them (a) and only fit three of them (b, c and d). It was updated to support four variables in 412d8d296c2f18ab0353d2b8d6c1cd3b6011ccaf, so the code can be modified to fit all four variables now.

[zifredder]

Is it possible to implement a variant of Fit.Curve that works with a function with an arbitrary number of parameters? It could have the following signature: public static double[] Curve (double[] x, double[] y, Func<double[], double, double> f, double[] initialGuesses, double tolerance, int maxIterations) With the following behavior (or so -- I hope you get the idea):

public static double foo(double x, double[] parameters) {
    // This is a model with 2 parameters
    return f(x, parameters[0]) + g(x, parameters[1]);
}
public static double bar(double x, double[] parameters) {
    // This is a model with 3 parameters
    return f(x, parameters[0]) + g(x, parameters[1]) + h(x, parameters[3]);
}

var optimizedParametersFoo = Fit.Curve(x, y, (p, x) => foo(x, p), initialGuessesFoo);
var optimizedParametersBar = Fit.Curve(x, y, (p, x) => bar(x, p), initialGuessesBar);

where initialGuessesFoo and optimizedParametersFoo both hold two values and initialGuessesBar and optimizedParametersBar both hold three values. I tried around myself, but unfortunately I didn't get far because I'm new to C# and have no clue how the modules of your library work together.