Curve Fitting Toolbox: Methods guide

Mathematical formulation of nonlinear least squares, Trust Region Reflective optimization, robust loss functions, covariance estimation, confidence bands, and residual diagnostics.

← Curve Fitting Toolbox overview · Getting started · Modeling guide · Examples · Reference

Mathematical formulation

The Curve Fitting Toolbox solves bounded, weighted, robust nonlinear regression problems of the form:

\[\min_{\mathbf{p} \in \mathbb{R}^p} \sum_{i=1}^n \rho\left(\frac{y_i - f(x_i, \mathbf{p})}{\sigma_{\text{eff}, i}}\right) \quad \text{subject to} \quad l_j \le p_j \le u_j \quad (j=1, \dots, p)\]

where: - \(\mathbf{p} = [p_1, p_2, \dots, p_p]^T\) is the vector of model parameters to be estimated. - \(f(x, \mathbf{p})\) is the nonlinear model response function. - \((x_i, y_i)\) are the measured independent and dependent observations (\(i=1, \dots, n\)). - \(\sigma_{\text{eff}, i} = \frac{\sigma_i}{\sqrt{w_i}}\) is the effective observation standard deviation combining measurement uncertainty \(\sigma_i\) and relative statistical weight \(w_i\). - \(\rho(z)\) is the loss penalty function (\(z = r_i / \sigma_{\text{eff}, i}\)). - \(l_j\) and \(u_j\) represent strict lower and upper box bounds on parameter \(p_j\).

When parameters are locked using Fixed = TRUE, the optimization vector is partitioned into free parameters \(\mathbf{p}_{\text{free}} \in \mathbb{R}^{p_{\text{free}}}\) and fixed constants \(\mathbf{p}_{\text{fixed}}\), optimizing only over \(\mathbf{p}_{\text{free}}\).

Optimization algorithms: TRF vs Levenberg-Marquardt

The toolbox utilizes SciPy’s implementation of the Trust Region Reflective (TRF) algorithm (method='trf' in scipy.optimize.curve_fit):

Why Trust Region Reflective (TRF)?

  1. Strict Box Constraints: Unlike classical Levenberg-Marquardt (which is inherently unconstrained), TRF handles explicit upper and lower parameter bounds (\(l_j \le p_j \le u_j\)) directly. When an iteration step crosses a parameter boundary, TRF uses an affine transformation of the coordinate space and reflects the search path back into the feasible interior.
  2. Robust Loss Compatibility: TRF seamlessly minimizes non-quadratic robust loss functions (Huber, Soft L1, Cauchy) by reformulating the objective as an iteratively reweighted least squares trust-region subproblem.
  3. Global Convergence: By maintaining a dynamic trust-region radius \(\Delta_k\), TRF guarantees convergence to a local minimum even when initial guesses are distant from the true solution or when the Jacobian is locally rank-deficient.

Comparison: TRF vs Levenberg-Marquardt (LM)

Feature Trust Region Reflective (TRF) Levenberg-Marquardt (LM)
Parameter Bounds Supports strict lower and upper bounds (\(l_j \le p_j \le u_j\)) Unconstrained only (\(-\infty < p_j < \infty\))
Robust Loss Functions Supports Linear, Huber, Soft L1, and Cauchy Restricted to standard quadratic loss (Linear)
Under-determined Systems Handles rank-deficient Jacobians via trust-region damping Requires artificial regularization or fails to invert
Computational Cost Solves trust-region subproblem per step; slightly higher per iteration Solves linear system per step; faster per iteration
Default in SciPy Automatic when bounds or non-linear losses are provided Restricted to unbounded unweighted least squares

Robust loss functions

In real-world experiments, measurement errors often follow heavy-tailed distributions rather than pure Gaussians. A single outlier can disproportionately pull standard least squares estimates away from true physical values. The toolbox provides four loss penalty functions \(\rho(z)\), where \(z = \frac{y_i - f(x_i, \mathbf{p})}{\sigma_{\text{eff}, i}}\):

Loss Function Comparison: rho(z) penalty growth
  Linear:  quadratic growth (z^2) -> extreme outlier sensitivity
  Huber:   quadratic for |z| <= 1, linear for |z| > 1 -> resistant to moderate spikes
  Soft L1: smooth hyperbolic curve 2*(sqrt(1+z^2) - 1) -> smooth transition
  Cauchy:  logarithmic growth ln(1 + z^2) -> strong suppression of severe anomalies

Mathematical definitions:

  1. Linear (Standard Least Squares): \[\rho(z) = z^2\] \[\psi(z) = \rho'(z) = 2z\] Influence grows without bound. Optimal when errors are strictly normally distributed (\(\mathcal{N}(0, \sigma^2)\)).

  2. Huber Loss: \[\rho(z) = \begin{cases} z^2 & \text{if } |z| \le 1 \\ 2|z| - 1 & \text{if } |z| > 1 \end{cases}\] \[\psi(z) = \begin{cases} 2z & \text{if } |z| \le 1 \\ 2\,\text{sgn}(z) & \text{if } |z| > 1 \end{cases}\] Limits the influence of residuals greater than \(1\sigma\) to a constant value, providing robust resistance against contamination up to 15–20% outliers.

  3. Soft L1 Loss: \[\rho(z) = 2\left(\sqrt{1 + z^2} - 1\right)\] \[\psi(z) = \frac{2z}{\sqrt{1 + z^2}}\] A smooth, infinitely differentiable approximation to absolute value (L1) loss. Residual influence saturates asymptotically at \(\pm 2\).

  4. Cauchy (Lorentzian) Loss: \[\rho(z) = \ln(1 + z^2)\] \[\psi(z) = \frac{2z}{1 + z^2}\] The influence function \(\psi(z)\) actually decreases toward zero for large residuals (\(|z| \to \infty\)). This provides strong redescending robustness, effectively ignoring catastrophic measurement spikes or instrumental failures.

Asymptotic covariance matrix and parameter inference

Upon finding the optimal parameter vector \(\mathbf{p}^*\), the toolbox computes statistical uncertainties via the linearized Jacobian matrix:

\[J_{ij} = -\left.\frac{\partial f(x_i, \mathbf{p})}{\partial p_j}\right|_{\mathbf{p}^*}\]

For free parameters, the unscaled covariance matrix is given by:

\[\mathbf{C}_{\text{unscaled}} = \left(\mathbf{J}_{\text{free}}^T \mathbf{W} \mathbf{J}_{\text{free}}\right)^{-1}\]

where \(\mathbf{W} = \text{diag}\left(\frac{w_1}{\sigma_1^2}, \dots, \frac{w_n}{\sigma_n^2}\right)\). When observational sigmas are relative rather than absolute (the standard setting), the covariance matrix is scaled by the estimated residual variance (reduced chi-squared):

\[s^2 = \frac{\text{SSE}}{n - p_{\text{free}}} = \frac{\sum_{i=1}^n (y_i - f(x_i, \mathbf{p}^*))^2}{n - p_{\text{free}}}\] \[\mathbf{C} = s^2 \cdot \mathbf{C}_{\text{unscaled}}\]

Parameter statistics

  • Standard Error: For free parameter \(j\), the standard error is the square root of the corresponding diagonal variance: \[\text{SE}(p_j) = \sqrt{C_{jj}}\] For fixed parameters, \(\text{SE}(p_k) \equiv 0.0\).
  • Student-t Statistic: \[t_j = \frac{\hat{p}_j}{\text{SE}(p_j)}\]
  • p-Value: \[p = 2 \cdot \left(1 - F_t(|t_j|; \nu)\right)\] where \(F_t\) is the cumulative distribution function of Student’s t-distribution with \(\nu = n - p_{\text{free}}\) residual degrees of freedom.
  • Asymptotic \((1 - \alpha)\) Confidence Interval: \[\left[\hat{p}_j - t_{\alpha/2, \nu} \cdot \text{SE}(p_j), \quad \hat{p}_j + t_{\alpha/2, \nu} \cdot \text{SE}(p_j)\right]\]

Pointwise curve confidence bands

In addition to individual parameter confidence intervals, the toolbox generates pointwise confidence bands along the fitted curve \(f(x, \hat{\mathbf{p}})\) across the evaluation grid.

Using the Delta method, the variance of the predicted response \(\hat{y}(x) = f(x, \hat{\mathbf{p}})\) is obtained by projecting the parameter covariance matrix through the local gradient vector:

\[\mathbf{g}(x) = \nabla_{\mathbf{p}} f(x, \hat{\mathbf{p}}) = \left[\frac{\partial f(x, \mathbf{p})}{\partial p_1}, \dots, \frac{\partial f(x, \mathbf{p})}{\partial p_p}\right]^T\] \[\text{Var}(\hat{y}(x)) = \mathbf{g}(x)^T \mathbf{C}_{\text{full}} \mathbf{g}(x)\] \[\text{SE}(\hat{y}(x)) = \sqrt{\text{Var}(\hat{y}(x))}\]

The pointwise \((1 - \alpha)\) confidence limits published on the Fitted Curve sheet are:

\[\text{CI}_{\text{lower}}(x) = \hat{y}(x) - t_{\alpha/2, \nu} \cdot \text{SE}(\hat{y}(x))\] \[\text{CI}_{\text{upper}}(x) = \hat{y}(x) + t_{\alpha/2, \nu} \cdot \text{SE}(\hat{y}(x))\]

In the task-pane plot, this is displayed as a shaded confidence envelope (\(p=0.95\)), illustrating where model predictions are most certain (typically near the centroid of measured data) and where prediction uncertainty expands (at domain boundaries or extrapolation regions).

Numerical derivatives

All derivatives required for Jacobian construction, Delta-method covariance projection, and grid curve slopes are evaluated numerically using central finite differences with step size \(h = 10^{-6}\):

1. Parameter gradient:

\[\frac{\partial f(x, \mathbf{p})}{\partial p_j} \approx \frac{f(x, \mathbf{p} + h \mathbf{e}_j) - f(x, \mathbf{p} - h \mathbf{e}_j)}{2h}\]

2. Fitted curve first derivatives:

Published in column Derivative of Fitted Curve: - Interior points (\(1 < i < N\)): Central difference: \[f'(x_i) \approx \frac{y_{i+1} - y_{i-1}}{x_{i+1} - x_{i-1}}\] - Left boundary (\(i=1\)): Forward difference: \[f'(x_1) \approx \frac{y_2 - y_1}{x_2 - x_1}\] - Right boundary (\(i=N\)): Backward difference: \[f'(x_N) \approx \frac{y_N - y_{N-1}}{x_N - x_{N-1}}\]

Diagnostics Suite metrics and thresholds

The 7-check Diagnostics Suite evaluates fit validity across seven rigorous statistical criteria:

1. Model validation

  • Formula: \(n \ge p_{\text{free}} + 1\)
  • Severity: PASS if \(n > p_{\text{free}}\), else FAIL.
  • Detail: Verifies that sufficient independent observations exist to avoid an under-determined system.

2. Convergence and boundary proximity

  • Formula: Check solver termination flag and boundary distance \(|p_j - l_j| \le 10^{-4}\) or \(|p_j - u_j| \le 10^{-4}\) for all \(j \in \text{free}\).
  • Severity: PASS if converged in interior; WARN if any free parameter hits its bound; FAIL if solver failed.
  • Detail: A parameter stuck at its bound indicates that the physical optimum lies outside the declared range or that the model is overparameterized.

3. Goodness of fit (\(R^2\))

  • Formula: \(R^2 = 1 - \frac{\text{SSE}}{\text{SST}}\) where \(\text{SST} = \sum_{i=1}^n (y_i - \bar{y})^2\).
  • Severity:
    • \(R^2 \ge 0.95\): PASS (high explanatory power)
    • \(0.80 \le R^2 < 0.95\): WARN (moderate explanatory power or high noise)
    • \(R^2 < 0.80\): FAIL (poor fit; model equation may be inappropriate)

4. Scaled Jacobian conditioning (\(\kappa\))

  • Formula: To avoid false ill-conditioning caused by disparate parameter units, the Jacobian columns are normalized: \[\mathbf{D} = \text{diag}(\|\mathbf{J}_{*, 1}\|_2, \dots, \|\mathbf{J}_{*, p}\|_2), \quad \mathbf{J}_{\text{scaled}} = \mathbf{J} \mathbf{D}^{-1}\] \[\kappa = \text{cond}(\mathbf{J}_{\text{scaled}}) = \frac{\sigma_{\max}(\mathbf{J}_{\text{scaled}})}{\sigma_{\min}(\mathbf{J}_{\text{scaled}})}\]
  • Severity:
    • \(\kappa < 1000\): PASS (well-conditioned; parameters are statistically independent)
    • \(1000 \le \kappa < 10000\): WARN (moderate collinearity; parameter estimates may have high correlation)
    • \(\kappa \ge 10000\): FAIL (severe ill-conditioning; model is overparameterized or data lacks dynamic range)

5. Shapiro-Wilk residual normality

  • Formula: \[W = \frac{\left(\sum_{i=1}^n a_i r_{(i)}\right)^2}{\sum_{i=1}^n (r_i - \bar{r})^2}\]
  • Severity: PASS if \(p \ge 0.05\); WARN if \(p < 0.05\).
  • Detail: Tests whether residuals conform to a normal distribution. Systematic non-normality suggests missing terms or heteroscedastic variance.

6. Durbin-Watson autocorrelation (\(d\))

  • Formula: \[d = \frac{\sum_{i=2}^n (r_i - r_{i-1})^2}{\sum_{i=1}^n r_i^2}\]
  • Severity:
    • \(1.5 \le d \le 2.5\): PASS (no evidence of serial correlation)
    • \(1.0 \le d < 1.5\) or \(2.5 < d \le 3.0\): WARN (mild positive or negative correlation)
    • \(d < 1.0\) or \(d > 3.0\): FAIL (severe autocorrelation; systematic curvature missing from model)

7. Standardized outlier detection

  • Formula: Standardized residual \(z_i = \frac{|r_i|}{\text{RMSE}}\).
  • Severity:
    • Auto mode: PASS if all \(z_i \le 2.5\); WARN if any \(z_i > 2.5\).
    • Strict mode: PASS if all \(z_i \le 2.0\); WARN if any \(z_i > 2.0\).
    • Disabled mode: PASS (detection bypassed).

Comparison with MATLAB Curve Fitting Toolbox

The Boardflare Curve Fitting Toolbox provides analytical capabilities analogous to the commercial MATLAB Curve Fitting Toolbox™ without requiring external licenses or MATLAB installations.

Command and API mapping

MATLAB Curve Fitting Toolbox Boardflare Curve Fitting Toolbox Functional equivalence
fit(x, y, 'rat11') / fittype Catalog michaelis_menten, exponential_decay, etc. Select parametric equation family
fitoptions('Method', 'NonlinearLeastSquares') Trust Region Reflective via scipy.optimize.curve_fit Bounded nonlinear least squares solver
opts.Lower = [...], opts.Upper = [...] Parameters!LowerBound, Parameters!UpperBound Strict parameter box constraints
opts.Exclude = ... / opts.Robust = 'LAR' Start!Loss function (Huber, Cauchy, Soft L1) Robust regression suppressing outliers
opts.Weights = w Data!Weight, Data!Sigma Pairwise observation weighting (\(\sigma_{\text{eff}} = \sigma / \sqrt{w}\))
confint(fitresult, 0.95) Results!CI_Lower, Results!CI_Upper Parameter 95% asymptotic confidence intervals
predint(fitresult, xgrid, 0.95, 'functional') Fitted Curve!CI_Lower, Fitted Curve!CI_Upper Pointwise 95% curve confidence bands
differentiate(fitresult, xgrid) Fitted Curve!Derivative Analytical / finite difference first derivatives
Residual plots & Durbin-Watson Diagnostics sheet 7-check automated suite Automated residual diagnostics and condition checks

Capability comparison

Dimension MATLAB Curve Fitting Toolbox Boardflare Curve Fitting Toolbox
Environment Standalone desktop application requiring MATLAB runtime Runs natively inside standard Microsoft Excel workbooks
Licensing Commercial per-seat MATLAB and Toolbox license required Open-source Python ecosystem (SciPy, NumPy, Pandas)
Data Flow Import/export from .mat files or workspace variables Live spreadsheet cells; updates when workbook cells change
Parameter Locking Requires reconfiguring fitoptions bounds (\(l_j = u_j\)) Dedicated Fixed = TRUE column with exact DoF preservation
Residual Diagnostics Interactive GUI inspection plots 7 automated PASS/WARN/FAIL checks published directly to Excel
Sharability Other users need MATLAB or compiled executables Share standard .xlsx with collaborators; viewable anywhere