Reineur Solver: solve real-world problems with confidence

Curve Fitting

Curve fitting is a mathematical technique used to construct a curve that best fits a series of data points. It is widely applied in data analysis, statistics, and machine learning to model relationships between variables.

Types of Curve Fitting:

Linear Regression: Fits a straight line to the data points.
Polynomial Regression: Fits a polynomial curve of degree n to the data points.
Nonlinear Regression: Fits a nonlinear model to the data points.

Example: Polynomial Curve Fitting

Given a set of data points, we can fit a polynomial curve using least squares optimization.

\min_{\mathbf{p}} \sum_{i=1}^{n} (y_i - P(x_i; \mathbf{p}))^2

Example 1C#

// Sample data points
double[] xData = [1, 2, 3, 4, 5];
double[] yData = [2.2, 3.0, 3.2, 2.5, 1.1];
// Fit a polynomial of degree 2
int degree = 2;
var coefficients = Polyfit(xData, yData, degree);
Scatter(xData, yData, "*", 15); HoldOn();
// Generate fitted curve
double[] xFit = Linspace(1, 5, 100);
double[] yFit = Polyval(coefficients, xFit);
Plot(xFit, yFit, Linewidth: 2);
SaveAs("Polynomial_Fitting.png");
CloseFig();

OutputFrom the book

Example: Fourier Series Fitting

Evaluating a Fourier series numerically involves transforming an infinite sum of trigonometric terms into a computationally stable, finite calculation while controlling truncation errors, floating-point precision loss, and spectral artifacts.

Mathematical Formulation:

A truncated Fourier series approximating a periodic function $f(x)$ on the interval $[-\pi, \pi]$ with $N$ harmonics is defined as:

f_N(x) = \frac{a_0}{2} + \sum_{n=1}^{N} \left( a_n \cos(nx) + b_n \sin(nx) \right)

In complex exponential form, which is computationally convenient for many numerical implementations, the series is expressed as:

f_N(x) = \sum_{n=-N}^{N} c_n e^{i n x}

where the complex coefficients $c_n$ relate to the real coefficients via:

c_0 = \frac{a_0}{2}, \quad c_n = \frac{a_n - i b_n}{2}, \quad c_{-n} = \frac{a_n + i b_n}{2}

Example 2C#

ColVec x = Linspace(-10, 10, 1001);
ColVec Rect = Sign(Sin(x));
Plot(x, Rect, Linewidth: 2); HoldOn();
var fourier = Plot(x, 0 * x, "r", Linewidth: 2);
Axis([x[0], x[^1], -1.5, 1.5]);
byte[] Animfun(int N)
{
    Matrix A = Zeros(1001, 2 * N + 3);
    A[.., 0] = Ones(1001);
    for (int n = 1; n <= (N + 1); n++)
    {
        A[.., 2 * n-1] = Cos(n * x); 
        A[.., 2 * n] = Sin(n * x);
    }
    ColVec p = Mldivide(A, Rect);
    fourier.Ydata = A * p;
    return GetFrame();
}
AnimationMaker(Animfun, "FourierFitting.gif", 5, 100);
CloseFig();

Example: Bi-Exponential Curve Fitting

This exercise covers non-linear parameter estimation using least-squares optimization to fit a bi-exponential model to noisy data while visualizing optimizer convergence.

The objective is to fit data points $(x_d, y_d)$ to a bi - exponential model:

f(x; \theta) = \theta_2 e^{\theta_0 x} + \theta_3 e^{\theta_1 x}

where math:\theta = [\theta_0, \theta_1, \theta_2, \theta_3]^T represents the unknown parameters.

Find $\hat{\theta}$ minimizing the sum of squared residuals:

\hat{\theta} = \arg\min_{\theta} \sum_{d=1}^D (y_d - f(x_d; \theta))^2

Example 3C#

ColVec noise, weight = new double[100]; double[] x0;
static ColVec fun(ColVec x, ColVec xdata) => x[2] * Exp(x[0] * xdata) + x[3] * Exp(x[1] * xdata);
ColVec xdata = Linspace(0, 1); noise = Rand(xdata.Numel);
ColVec ydata = fun(x0 = [-4, -5, 4, -4], xdata) + 0.02 * noise;
x0 = [-1, -2, 1, -1]; weight[xdata < 0.5] = 1;
var opts = OptimSet(Display: true, MaxIter: 200, StepTol: 1e-6, OptimalityTol: 1e-6);
var ans = Lsqcurvefit(fun, x0, xdata, ydata, options: opts);
AnimateHistory(fun, xdata, ydata, ans.history, "Bi_Exponential_Fitting.gif");
CloseFig();

OutputFrom the book

                                           Norm of      First-order 
Iteration   Func-count       Resnorm          step       optimality 
    0            5          1.7065e0                       3.5565e0 
    1           11         7.8717e-1     2.5350e-1         1.4636e0 
    2           17         5.1948e-1     2.3489e-1        2.9532e-1 
    3           23         4.4230e-1     5.2110e-1        1.4316e-1 
    4           29         2.6309e-1      2.0028e0        6.8482e-1 
    5           35         1.0851e-1      2.9815e0        7.2106e-1 
    6           41         5.4577e-3     4.0791e-1        9.0544e-2 
    7           47         4.3846e-3     4.1729e-1        1.1143e-1 
    8           53         3.9921e-3     3.5063e-1        9.4193e-2 
    9           59         3.5093e-3     2.3213e-1        4.6637e-2 
   10           65         3.3666e-3     1.0678e-1        1.0289e-2 
   11           71         3.3594e-3     2.4916e-2        5.6414e-4 
   12           77         3.3594e-3     2.0379e-3        3.7770e-6 
   13           83         3.3594e-3     3.8937e-5        1.5011e-9

Example 4C#

ColVec xdata, ydata, times, y_est, filltime, sgy, filly, lower, upper;
double[] x_dat = [0.9, 1.5, 13.8, 19.8, 24.1, 28.2, 35.2, 60.3, 74.6, 81.3];
double[] y_dat = [455.2, 428.6, 124.1, 67.3, 43.2, 28.1, 13.1, -0.4, -1.3, -1.5];
xdata = x_dat; ydata = y_dat; times = Linspace(x_dat[0], x_dat[9]);
double[] x0 = [100, -1];
static ColVec fun(ColVec x, ColVec xdata) => x[0] * Exp(x[1] * xdata);
var opts = OptimSet(Display: true, MaxIter: 200, StepTol: 1e-6, OptimalityTol: 1e-6);
var ans = Lsqcurvefit(fun, x0, xdata, ydata, options: opts);
Scatter(xdata, ydata); HoldOn();
Plot(times, y_est = fun(ans.x, times), "r", Linewidth: 2);
filltime = Vcart(times, times.Reverse().ToList());
sgy = Interp1(xdata, ans.sigma_y, times);
lower = y_est - 20 * sgy; upper = y_est + 20 * sgy;
filly = Vcart(lower, upper.Reverse().ToList());
Fill(filltime, filly, "g", 0.2); HoldOff();
Axis([xdata.Min()-0.01*xdata.Range(), xdata.Max()+0.01*xdata.Range(),
ydata.Min()-0.1*ydata.Range(), ydata.Max()+0.1*ydata.Range()]);
SaveAs("CurveFitting.png");
AnimateHistory(fun, xdata, ydata, ans.history, "CurveFitting.gif");
CloseFig();

ColVec xdata, ydata, times, y_est, filltime, sgy, filly, lower, upper;

double[] x_dat = [0.9, 1.5, 13.8, 19.8, 24.1, 28.2, 35.2, 60.3, 74.6, 81.3];
double[] y_dat = [455.2, 428.6, 124.1, 67.3, 43.2, 28.1, 13.1, -0.4, -1.3, -1.5];
xdata = x_dat; ydata = y_dat; times = Linspace(x_dat[0], x_dat[9]);
double[] x0 = [100, -1];

static ColVec fun(ColVec x, ColVec xdata) => x[0] * Exp(x[1] * xdata);
var opts = OptimSet(Display: true, MaxIter: 200, StepTol: 1e-6, OptimalityTol: 1e-6);
var ans = Lsqcurvefit(fun, x0, xdata, ydata, options: opts);

Scatter(xdata, ydata); HoldOn();
Plot(times, y_est = fun(ans.x, times), "r", Linewidth: 2);
filltime = Vcart(times, times.Reverse().ToList());
sgy = Interp1(xdata, ans.sigma_y, times);
lower = y_est - 20 * sgy; upper = y_est + 20 * sgy;
filly = Vcart(lower, upper.Reverse().ToList());
Fill(filltime, filly, "g", 0.2); HoldOff();
Axis([xdata.Min()-0.01*xdata.Range(), xdata.Max()+0.01*xdata.Range(),
ydata.Min()-0.1*ydata.Range(), ydata.Max()+0.1*ydata.Range()]);
SaveAs("CurveFitting.png");
AnimateHistory(fun, xdata, ydata, ans.history, "CurveFitting.gif");
CloseFig();

OutputFrom the book

                                           Norm of      First-order 
Iteration   Func-count       Resnorm          step       optimality 
    0            3          3.5968e5                       2.8768e4 
    1            7          2.9148e5      4.5301e1         6.3631e4 
    2           11          1.4328e5      7.0536e1         1.8724e5 
    3           15          5.8838e4      8.1015e1         1.7583e5 
    4           19          2.1604e4      7.9171e1         1.3573e5 
    5           23          2.4371e3      8.1537e1         4.6492e4 
    6           27          6.2429e1      3.5477e1         8.8212e3 
    7           31          9.6405e0      5.5200e0         5.2344e2 
    8           35          9.5049e0     2.7383e-1         4.5771e0 
    9           39          9.5049e0     3.5902e-3        1.3319e-2 
   10           43          9.5049e0     9.0844e-6        5.6776e-6

Lsqcurvefit allows the use of constraints. 1. Seed data for reproducability

Example 5C#

int seed = 23;
Random rng = new(seed);
ColVec xdata, ydata, noise = Randn(100);
double[] xstar = [2, 4, 5, 0.5];
ColVec model(ColVec x, ColVec xdata) => x[0] + x[1] * Atan(xdata - x[2]) + x[3] * xdata;
xdata = Linspace(2,7); ydata = model(xstar, xdata) + noise/10;
Scatter(xdata, ydata, "ro");
Xlabel("x"); Ylabel("y"); 
SaveAs("Seeded_Curve_Fitting_Data.png");
CloseFig();

OutputFrom the book

Fitting with Linear constraint

Example 6C#

int seed = 23;
Random rng = new(seed);
ColVec xdata, ydata, noise = Randn(100);
double[] xstar = [2, 4, 5, 0.5], startpt = [1, 2, 3, 1];
ColVec model(ColVec x, ColVec xdata) => x[0] + x[1] * Atan(xdata - x[2]) + x[3] * xdata;
RowVec A = new double[] { -1, -1, 1, 1 };
ColVec fineq(ColVec x) => A * x; ColVec lb = Zeros(4), ub = 7 + lb;
xdata = Linspace(2, 7); ydata = model(xstar, xdata) + noise / 10;
var opts = OptimSet(Display: true, MaxIter: 200, StepTol: 1e-6, OptimalityTol: 1e-6);
var ans = Lsqcurvefit(model, startpt, xdata, ydata, fineq, null, lb, ub, options: opts);
Console.WriteLine($"x = {ans.x.T}");
Console.WriteLine($"c = {fineq(ans.x)}");
Scatter(xdata, ydata, "ro"); HoldOn();
Plot(xdata, ans.y_hat, "-b", Linewidth: 2);
Xlabel("x"); Ylabel("y");
Legend(["Measured Data", "Model Estimate"], UpperRight);
SaveAs("Example_of_CurveFitting_using_Lsqcurvefit_with_Linear_Inequality_Constraints.png");
CloseFig();

OutputFrom the book

                                            Norm of      First-order 
 Iteration   Func-count       Resnorm          step       optimality 
     0            5          1.5833e3                       1.5193e3 
     1           11          1.3013e3     2.1063e-1         1.3194e3 
     2           17          7.9592e2     4.8507e-1         8.9285e2 
     3           23          3.0032e2     8.0730e-1         3.3850e2 
     4           29          6.0406e1     9.6617e-1         6.9159e1 
     5           35          1.0621e1     6.5279e-1         3.4324e1 
     6           41          5.7907e0     3.8978e-1         2.2699e0 
     7           47          2.8720e0     5.4306e-1         3.3866e0 
     8           53          2.5536e0     5.6541e-1         1.8019e0 
     9           59          1.0647e0     2.7129e-1        3.2612e-1 
    10           65          1.0626e0     3.9289e-3        7.3238e-2 
    11           71          1.0532e0     1.5675e-1        5.5043e-2 
    12           77          1.0530e0     2.0308e-2        1.2136e-3 
    13           83          1.0530e0     8.7077e-4        2.2264e-5 
    14           89          1.0530e0     1.2797e-5        4.4444e-7 
x = 
   1.8356    3.9412    4.9886    0.5266

c =   -0.2616

Example_of_CurveFitting_using_Lsqcurvefit_with_Linear_Inequality_Constraints.png

Fitting with nonlinear constraint.

Example 7C#

int seed = 23;
Random rng = new(seed);
ColVec xdata, ydata, noise = Randn(100);
double[] xstar = [2, 4, 5, 0.5], startpt = [1, 2, 3, 1];
ColVec model(ColVec x, ColVec xdata) => x[0] + x[1] * Atan(xdata - x[2]) + x[3] * xdata;
ColVec fineq(ColVec x) => x[0] * x[0] + x[1] * x[1] - 16; ColVec lb = Zeros(4), ub = 7 + lb;
xdata = Linspace(2, 7); ydata = model(xstar, xdata) + noise / 10;
var opts = OptimSet(Display: true, MaxIter: 200, StepTol: 1e-6, OptimalityTol: 1e-6);
var ans = Lsqcurvefit(model, startpt, xdata, ydata, fineq, null, lb, ub, options: opts);
Console.WriteLine($"x = {ans.x.T}");
Console.WriteLine($"c = {fineq(ans.x)}");
Scatter(xdata, ydata, "ro"); HoldOn();
Plot(xdata, ans.y_hat, "-b", Linewidth: 2);
Xlabel("x"); Ylabel("y");
Legend(["Measured Data", "Model Estimate"], UpperRight);
SaveAs("Example_of_CurveFitting_using_Lsqcurvefit_with_NonLinear_Inequality_Constraints.png");
CloseFig();

OutputFrom the book

                                            Norm of      First-order 
 Iteration   Func-count       Resnorm          step       optimality 
     0            5          1.5708e3                       1.5056e3 
     1           11          1.2895e3     2.1199e-1         1.3055e3 
     2           17          7.8528e2     4.8911e-1         8.7877e2 
     3           23          2.8978e2     8.1787e-1         3.2570e2 
     4           29          5.3728e1     9.6956e-1         6.9713e1 
     5           35          5.0915e0     6.2595e-1         3.8177e1 
     6           41          2.3765e0     1.9906e-1         4.3527e0 
     7           47          2.1520e0     1.6511e-1         1.0622e0 
     8           53          1.7709e0     3.6655e-1         1.3801e0 
     9           59          1.2994e0     6.5160e-1         1.4022e0 
    10           68          1.2686e0     6.0764e-2        6.0069e-1 
    11           75          1.2424e0     5.6062e-2        5.6590e-1 
    12           83          1.2347e0     1.7305e-2        5.4868e-1 
    13           91          1.2324e0     5.4304e-3        5.4410e-1 
    14           99          1.2316e0     1.7131e-3        5.4274e-1 
    15          107          1.2314e0     5.4130e-4        5.4232e-1 
    16          114          1.2312e0     5.4085e-4        5.4193e-1 
    17          120          1.2309e0     3.5138e-4        3.6876e-1 
    18          126          1.2305e0     8.6562e-4        2.6907e-1 
    19          132          1.2300e0     1.5182e-3        2.3884e-1 
    20          138          1.2297e0     1.9533e-3        2.8957e-1 
    21          144          1.2296e0     1.2517e-3        3.3709e-1 
    22          150          1.2296e0     2.9346e-4        3.4901e-1 
    23          156          1.2296e0     2.2291e-5        3.4993e-1 
    24          162          1.2296e0     5.0231e-7        3.4995e-1 
x = 
   1.3386    3.7694    5.0124    0.6429

c =    0.0005

Example_of_CurveFitting_using_Lsqcurvefit_with_NonLinear_Inequality_Constraints.png

Takes a free account.