implementation_section.tex

% ============================================================
% Implementation section
% Covers both MATLAB and Python implementations of LinearRegression.
% Intended as a single \section{} or \subsection{} of a larger document.
%
% Required packages (add to preamble if not already present):
%   \usepackage{booktabs, array, listings, tikz, xcolor, textcomp}
%   \usetikzlibrary{shapes.geometric, arrows.meta}
%
% Required macro (add to preamble):
%   % Renders a quoted code string with upright single quotes, safe in tables.
%   % Usage: \mstr{}  ->  ''        \mstr{x}  ->  'x'        \mstr{WLS}  ->  'WLS'
%   \newcommand{\mstr}[1]{\texttt{\textquotesingle#1\textquotesingle}}
%
% Listing styles required:
%   \lstdefinestyle{matlab}{...}   % already defined in parent document
%   \lstdefinestyle{python}{...}   % defined below — add to preamble
%
% Suggested preamble addition for Python style:
%   \lstdefinestyle{python}{
%     language=Python, basicstyle=\ttfamily\small,
%     keywordstyle=\color{blue!70!black}\bfseries,
%     stringstyle=\color{orange!80!black},
%     commentstyle=\color{gray}\itshape,
%     showstringspaces=false, breaklines=true,
%     frame=single, numbers=left, numberstyle=\tiny\color{gray},
%     morekeywords={True,False,None}
%   }
% ============================================================

This section documents the implementations of the WLS and TLS estimators
derived in the preceding sections: a MATLAB implementation in
\texttt{LinearRegression.m} and a Python implementation in
\texttt{linear\_regression.py}.  Both build the design matrix automatically
from a symbolic model string and compute the complete uncertainty budget
described in Section~\ref{sec:Uncertainty}.

The Python version began as a direct port of the MATLAB one, and the two
produce the same numbers: same estimator, same covariance, same
$\max(1,\tilde\chi^2)$ inflation rule.  They differ in one respect only.
MATLAB draws the diagnostic figure from inside the fit call, whereas in Python
the plot is a separate public function and the plotting options have moved with
it --- so the two option lists are not interchangeable.

\noindent The following subsections document the shared design first, then the
language-specific calling conventions and output formats separately.

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Code Architecture}
% ─────────────────────────────────────────────────────────────────────────────

Each implementation consists of one public entry-point and a set of
subfunctions, summarised in Table~\ref{tab:subfunctions}.  The plotting
routine is where the two differ: in MATLAB it is a private subfunction called
from inside the fit, whereas in Python it is a public module-level function
called separately, so a figure can be regenerated from a stored result without
re-running the fit.

\begin{table}[h]
\centering
\caption{Functions in each implementation.  Python private names are prefixed
         with an underscore by convention.}
\label{tab:subfunctions}
\begin{tabular}{p{4.4cm}llll}
\toprule
Role & MATLAB name & Python name & MATLAB & Python \\
\midrule
Argument parsing, NaN removal, solver, uncertainty budget
  & \texttt{LinearRegression}     & \texttt{linear\_regression}      & public  & public \\
Symbolic differentiation $\to \mathbf{X}$, $\mathbf{U}$
  & \texttt{build\_design\_matrix}& \texttt{\_build\_design\_matrix} & private & private \\
Formatted console output
  & \texttt{print\_summary}       & \texttt{\_print\_summary}        & private & private \\
Two-panel diagnostic figure
  & \texttt{make\_plot}           & \texttt{make\_plot}              & private & public \\
\bottomrule
\end{tabular}
\end{table}

\noindent The top-level execution flow is shown below; it is identical
in both implementations.

\begin{center}
\begin{tikzpicture}[
  box/.style      = {rectangle, draw, rounded corners=3pt,
                     text width=72mm, minimum height=10mm,
                     align=center, font=\small},
  decision/.style = {diamond, draw, aspect=3.2,
                     text width=30mm, minimum height=10mm,
                     align=center, font=\small},
  arr/.style      = {-Stealth, thick}
]
\node[box]      at (0,  0.0) (in)    {Inputs: $x,\,y,\,u_y$\,[opt: $u_x$], modelStr, params, options};
\node[decision] at (0, -1.6) (det)   {4th arg a string?};
\node[box]      at (0, -3.2) (nan)   {Remove NaN rows};
\node[box]      at (0, -4.6) (dm)    {\texttt{build\_design\_matrix} $\;\to\; \mathbf{X},\,\mathbf{U}$};
\node[box]      at (0, -6.0) (solve) {Solve: WLS (closed-form) or TLS (iterative IRLS)};
\node[box]      at (0, -7.4) (unc)   {Uncertainty budget: $\mathrm{Cov},\;\tilde\chi^2,\;
                                       w_\mathrm{rms},\;u_\mathrm{total}^{(A/B)},\;\bar{u}$};
\node[box]      at (0, -8.8) (out)   {Pack output (struct / dataclass)};
\node[box]      at (0,-10.2) (disp)  {[opt]\;\texttt{print\_summary}\quad(MATLAB also: [opt]\;\texttt{make\_plot})};

\draw[arr] (in.south)    -- (det.north);
\draw[arr] (det.south)   -- node[right, font=\small, xshift=1mm]{yes $\to$ WLS} (nan.north);
\draw[arr] (det.east)    -- ++(20mm,0)
           |- node[near start, right, font=\small, xshift=1mm]{no $\to$ TLS} (nan.east);
\draw[arr] (nan.south)   -- (dm.north);
\draw[arr] (dm.south)    -- (solve.north);
\draw[arr] (solve.south) -- (unc.north);
\draw[arr] (unc.south)   -- (out.north);
\draw[arr] (out.south)   -- (disp.north);
\end{tikzpicture}
\end{center}

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Design Matrix Construction}
\label{sec:impl_dm}
% ─────────────────────────────────────────────────────────────────────────────

Both implementations construct $\mathbf{X}$ and $\mathbf{U}$ analytically
from the model string using symbolic differentiation.  MATLAB uses the
Symbolic Math Toolbox (\texttt{syms}, \texttt{diff}, \texttt{matlabFunction});
Python uses SymPy (\texttt{sympy.Symbol}, \texttt{sympy.diff},
\texttt{sympy.lambdify}).  For each parameter $\beta_k$:
\begin{align}
  X_{ik} &= \left.\frac{\partial f}{\partial \beta_k}\right|_{x_i}
             \quad \text{(symbolic diff, evaluated numerically)},\\[4pt]
  U_{ik} &= \left|\frac{\partial^2 f}{\partial \beta_k\,\partial x}
             \right|_{x_i} u(x_i)
             \quad \text{(GUM first-order propagation, Eq.~\ref{eq:sigmaEffUncorr})}.
\end{align}
No numerical differencing is used; both derivatives are exact.

Model strings use standard mathematical notation with \texttt{x} as the
independent variable.  MATLAB uses \texttt{\^{}} for exponentiation;
Python accepts both \texttt{\^{}} and \texttt{**} (the \texttt{\^{}}
character is silently converted to \texttt{**} before parsing).  The same
model string therefore works in both languages:

\begin{lstlisting}[style=matlab, caption={MATLAB model string.}]
modelStr = 'A + B*log(x) + D*(log(x))^3';
params   = ["A", "B", "D"];
\end{lstlisting}

\begin{lstlisting}[style=python, caption={Python model string (identical syntax accepted).}]
model_str = 'A + B*log(x) + D*(log(x))^3'  # ^ auto-converted to **
params    = ['A', 'B', 'D']
\end{lstlisting}

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Iterative Solver (TLS)}
% ─────────────────────────────────────────────────────────────────────────────

The TLS solver implements iteratively reweighted least squares (IRLS).
The solver is selected automatically when the fourth argument is numeric
(MATLAB) or is not a string (Python).  At each iteration the effective
variance (Eq.~\ref{eq:sigmaEffUncorr}) is recomputed from the current
$\hat{\vec\beta}$, the weight matrix updated, and a standard WLS solve
performed.  Convergence is declared when
\begin{equation}
  \|\hat{\vec\beta}^{(k+1)} - \hat{\vec\beta}^{(k)}\|
  < \varepsilon_{\mathrm{tol}}\bigl(1 + \|\hat{\vec\beta}^{(k)}\|\bigr),
  \qquad \varepsilon_{\mathrm{tol}} = 10^{-10}\ \text{(default)}.
\end{equation}
If the maximum iteration count is reached without convergence, a warning is
issued and the final iterate is returned.  The initial guess defaults to the
ordinary least-squares solution (no weighting); a custom starting point may
be supplied via the \texttt{beta0} option.

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Uncertainty Budget Implementation}
% ─────────────────────────────────────────────────────────────────────────────

Following the derivation in Section~\ref{sec:Chi2Inflation}, both
scatter-handling options are computed and stored side-by-side in the output,
allowing the user to select the appropriate one post-hoc:

\begin{center}
\begin{tabular}{lll}
\toprule
Field & Formula & When to use \\
\midrule
\texttt{uTotal\_chi2scaled} & $\sqrt{\mathbf{x}_i^\mathsf{T}\mathrm{Cov}_\mathrm{inf}\,\mathbf{x}_i}$ &
  $\tilde\chi^2 \gg 1$, scatter from underestimated $u$ \\
\texttt{uTotal\_wrmsAdded}  & $\sqrt{u_{\hat{y}_i}^2 + w_\mathrm{rms}^2}$ &
  Genuine independent noise floor \\
\texttt{uBar\_chi2scaled}   & Eq.~(\ref{eq:uCalMean}), Option~A & Single reportable number \\
\texttt{uBar\_wrmsAdded}    & Eq.~(\ref{eq:uCalMean}), Option~B & Single reportable number \\
\bottomrule
\end{tabular}
\end{center}

\noindent These field names, and their contents, are identical in both
implementations.  In particular \texttt{covBeta\_inf} follows the
inflate-only rule of Section~\ref{sec:Chi2Inflation} in both:
\begin{equation}
  \mathrm{Cov}_\mathrm{inf} = \max\!\left(1,\tilde\chi^2\right)\mathrm{Cov},
\end{equation}
so the bands widen when the residuals exceed the uncertainty model and hold at
$\mathrm{Cov}$ when they do not.  They are never shrunk below what the input
uncertainties support.

\medskip
\noindent In MATLAB, do \emph{not} pass \texttt{covBeta\_inf} to the plot while
also setting \texttt{includeScatter} to \texttt{true}, as this double-counts
the scatter.  Python's plot takes the result object and selects the bands
itself, so the conflict cannot arise.

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Display-Unit Conversion}
% ─────────────────────────────────────────────────────────────────────────────

An optional transform function $g(v)$ maps fit-space quantities to physical
display units before printing and plotting.  The local Jacobian is evaluated
by a relative-step finite difference,
\begin{equation}
  g'(v_i) \approx \frac{g\!\left(v_i(1+\varepsilon)\right) - g(v_i)}{v_i\,\varepsilon},
  \qquad \varepsilon = 10^{-6},
\end{equation}
giving $u_\mathrm{display} = |g'(v_i)|\,u_\mathrm{fitspace}$.
The relative step ensures good numerical conditioning regardless of the
magnitude of $v_i$.

% ─────────────────────────────────────────────────────────────────────────────
\subsection{MATLAB Interface}
% ─────────────────────────────────────────────────────────────────────────────

\subsubsection*{Dependencies}

Requires MATLAB with the Symbolic Math Toolbox.

\subsubsection*{Calling Convention}

The solver is selected by inspecting the type of the fourth argument.
If it is a string, the call is WLS; if it is numeric, the call is TLS.

\begin{lstlisting}[style=matlab]
% WLS (no x-uncertainty)
result = LinearRegression(x, y, uY, modelStr, params, Name, Value, ...)

% TLS (with x-uncertainty)
result = LinearRegression(x, y, uX, uY, modelStr, params, Name, Value, ...)
\end{lstlisting}

\subsubsection*{Name-Value Options}

\begin{tabular}{llp{6.2cm}}
\toprule
Name & Default & Description \\
\midrule
\texttt{tol}            & \texttt{1e-10}  & TLS convergence tolerance \\
\texttt{maxIter}        & \texttt{1000}   & TLS maximum iterations \\
\texttt{beta0}          & OLS             & TLS initial guess ($p\times 1$) \\
\texttt{printSummary}   & \texttt{true}   & Print coefficient table and $\bar{u}$ to console \\
\texttt{makePlot}       & \texttt{true}   & Generate two-panel diagnostic figure \\
\texttt{transformFcn}   & \texttt{@(v)v}  & Fit-space $\to$ display-unit mapping \\
\texttt{displayUnit}    & \mstr{}     & Unit label for display axis \\
\texttt{residScale}     & \texttt{1}      & Residual multiplier (e.g.\ 1000 for milli-units) \\
\texttt{residUnit}      & \mstr{}     & Residual axis unit label \\
\texttt{xLabel}         & \mstr{x}    & Top-panel $x$-axis label \\
\texttt{yLabel}         & \mstr{y}    & Top-panel $y$-axis label \\
\texttt{xBotLabel}      & \texttt{yLabel} & Bottom-panel $x$-axis label \\
\texttt{plotTitle}      & \mstr{}     & Figure title \\
\texttt{nGrid}          & \texttt{200}    & Points in dense fit curve \\
\texttt{xTopData}       & $x$             & Override top-panel $x$-axis data \\
\texttt{xBotData}       & $y_\text{disp}$ & Override bottom-panel $x$-axis data \\
\texttt{includeScatter} & \texttt{true}   & Add $w_\mathrm{rms}$ to CI bands (Option~B);
                                            set \texttt{false} for Option~A \\
\bottomrule
\end{tabular}

\subsubsection*{Output Struct}

\begin{tabular}{lll}
\toprule
Field & Size & Description \\
\midrule
\texttt{beta}              & $p\times 1$ & $\hat{\vec\beta}$ \\
\texttt{uBeta}             & $p\times 1$ & $u(\hat\beta_j)$ from \texttt{covBeta} \\
\texttt{uBeta\_inf}        & $p\times 1$ & $u(\hat\beta_j)$ from \texttt{covBeta\_inf} \\
\texttt{covBeta}           & $p\times p$ & Covariance, input-$u$ based (Eq.~\ref{eq:CovTLSgeneral}) \\
\texttt{covBeta\_inf}      & $p\times p$ & $\tilde\chi^2$-inflated covariance \\
\texttt{chi2r}             & scalar      & $\tilde\chi^2$ \\
\texttt{wrms}              & scalar      & $w_\mathrm{rms}$ (fit-space units) \\
\texttt{uFit}              & $n\times 1$ & $u_{\hat y_i}$ from \texttt{covBeta} (Eq.~\ref{eq:uFit}) \\
\texttt{uFit\_inf}         & $n\times 1$ & $u_{\hat y_i}$ from \texttt{covBeta\_inf} \\
\texttt{uTotal\_chi2scaled}& $n\times 1$ & Per-point $u$, Option~A \\
\texttt{uTotal\_wrmsAdded} & $n\times 1$ & Per-point $u$, Option~B \\
\texttt{uBar\_chi2scaled}  & scalar      & $\bar{u}$ summary, Option~A (fit-space) \\
\texttt{uBar\_wrmsAdded}   & scalar      & $\bar{u}$ summary, Option~B (fit-space) \\
\texttt{uEff}              & $n\times 1$ & $u_{\mathrm{eff},i}$ \\
\texttt{X}, \texttt{uX}   & $n\times p$ & Design matrix and its uncertainty \\
\texttt{Xhat}              & $n\times p$ & Corrected design matrix (TLS; $=\mathbf{X}$ for WLS) \\
\texttt{x\_valid}, \texttt{y\_valid} & $n\times 1$ & Data after NaN removal \\
\texttt{method}            & string      & \mstr{WLS} or \mstr{TLS} \\
\texttt{modelStr}          & string      & Model string as passed \\
\texttt{params}            & string array& Parameter names as passed \\
\texttt{n}, \texttt{p}    & scalar      & Observations used; parameters \\
\bottomrule
\end{tabular}

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Python Interface}
% ─────────────────────────────────────────────────────────────────────────────

\subsubsection*{Dependencies}

\begin{lstlisting}[style=python]
pip install numpy sympy          # required
pip install matplotlib           # optional: only needed for make_plot
\end{lstlisting}

\noindent \texttt{matplotlib} is imported lazily inside \texttt{make\_plot},
which also forces the \texttt{Agg} backend so no window is ever opened.  The
module therefore imports cleanly without it, and \texttt{make\_plot} returns
\texttt{None} rather than raising if it is missing.

\subsubsection*{Calling Convention}

The dispatch logic is identical to MATLAB: the fourth argument is inspected
at runtime and the solver is chosen accordingly.  Unlike MATLAB, the fit does
\emph{not} draw a figure --- call \texttt{make\_plot} separately.

\begin{lstlisting}[style=python]
from linear_regression import linear_regression, make_plot

# WLS (no x-uncertainty)
result = linear_regression(x, y, uY, model_str, params, **kwargs)

# TLS (with x-uncertainty)
result = linear_regression(x, y, uX, uY, model_str, params, **kwargs)

# Plot, whenever you want, from a stored result
fig = make_plot(result, resid_scale=1000, resid_unit='mK')
fig.savefig('fit.png', dpi=150)   # caller saves and closes
\end{lstlisting}

All options are passed as Python keyword arguments rather than MATLAB
name-value pairs.  Option names use \texttt{snake\_case} throughout.

\subsubsection*{Keyword Arguments to \texttt{linear\_regression}}

\begin{tabular}{llp{6.2cm}}
\toprule
Keyword & Default & Description \\
\midrule
\texttt{tol}             & \texttt{1e-10}  & TLS convergence tolerance \\
\texttt{max\_iter}       & \texttt{1000}   & TLS maximum iterations \\
\texttt{beta0}           & OLS             & TLS initial guess (\texttt{ndarray}, $p\times 1$) \\
\texttt{print\_summary}  & \texttt{True}   & Print coefficient table and $\bar{u}$ to console \\
\texttt{transform\_fcn}  & \texttt{None}   & Fit-space $\to$ display-unit mapping (callable) \\
\texttt{display\_unit}   & \mstr{}     & Unit label for display axis \\
\texttt{resid\_scale}    & \texttt{1.0}    & Residual multiplier (e.g.\ 1000 for milli-units) \\
\texttt{resid\_unit}     & \mstr{}     & Residual axis unit label \\
\texttt{x\_label}        & \mstr{x}    & Label carried into the printed summary \\
\texttt{y\_label}        & \mstr{y}    & Label carried into the printed summary \\
\texttt{plot\_title}     & \mstr{}     & Title carried into the printed summary \\
\bottomrule
\end{tabular}

\subsubsection*{Keyword Arguments to \texttt{make\_plot}}

All arguments after \texttt{result} are keyword-only.  The function returns the
\texttt{matplotlib} \texttt{Figure}, or \texttt{None} if \texttt{matplotlib} is
not installed; saving and closing it is the caller's responsibility.

\begin{tabular}{llp{6.2cm}}
\toprule
Keyword & Default & Description \\
\midrule
\texttt{transform\_fcn}  & \texttt{None}   & Fit-space $\to$ display-unit mapping (callable) \\
\texttt{resid\_scale}    & \texttt{1.0}    & Residual multiplier (e.g.\ 1000 for milli-units) \\
\texttt{resid\_unit}     & \mstr{}     & Residual axis unit label \\
\texttt{x\_label}        & \mstr{x}    & Top-panel $x$-axis label \\
\texttt{y\_label}        & \mstr{y}    & Top-panel $y$-axis label \\
\texttt{x\_bot\_label}   & \mstr{}     & Bottom-panel $x$-axis label (defaults to \texttt{y\_label}) \\
\texttt{plot\_title}     & \mstr{}     & Figure title \\
\texttt{n\_grid}         & \texttt{200}    & Points in dense fit curve \\
\texttt{x\_top\_data}    & \texttt{None}   & Override top-panel $x$-axis data \\
\texttt{x\_bot\_data}    & \texttt{None}   & Override bottom-panel $x$-axis data \\
\texttt{k}               & \texttt{2.0}    & Coverage factor for the plotted bands and error bars \\
\texttt{show\_inflated}  & \texttt{True}   & Also draw the $\tilde\chi^2$-inflated band over the
                                             non-inflated one; \texttt{False} draws only the latter \\
\texttt{floor\_disp}     & \texttt{0.0}    & Systematic floor added to the displayed bands \\
\bottomrule
\end{tabular}

\medskip
\noindent Both bands are drawn at the same coverage factor $k$ and overlap
semi-transparently, so the inflated and non-inflated intervals are visible at
once.  Note that $k=1$ is deliberately never used by this figure; the bands
shown in Figure~\ref{fig:ThermistorFit} were produced by the MATLAB routine,
which plots at $k=1$.

\subsubsection*{Output Dataclass}

The function returns a \texttt{RegressionResult} dataclass instance.
All fields support dot-access (\texttt{result.beta}, \texttt{result.chi2r},
etc.), matching MATLAB struct syntax.  Field names are identical to the
MATLAB struct with one exception: \texttt{modelStr} $\to$ \texttt{model\_str}.

\begin{tabular}{lll}
\toprule
Field & Size & Description \\
\midrule
\texttt{beta}              & $p\times 1$ & $\hat{\vec\beta}$ \\
\texttt{uBeta}             & $p\times 1$ & $u(\hat\beta_j)$ from \texttt{covBeta} \\
\texttt{uBeta\_inf}        & $p\times 1$ & $u(\hat\beta_j)$ from \texttt{covBeta\_inf} \\
\texttt{covBeta}           & $p\times p$ & Covariance, input-$u$ based (Eq.~\ref{eq:CovTLSgeneral}) \\
\texttt{covBeta\_inf}      & $p\times p$ & $\tilde\chi^2$-inflated covariance \\
\texttt{chi2r}             & scalar      & $\tilde\chi^2$ \\
\texttt{wrms}              & scalar      & $w_\mathrm{rms}$ (fit-space units) \\
\texttt{uFit}              & $n\times 1$ & $u_{\hat y_i}$ from \texttt{covBeta} (Eq.~\ref{eq:uFit}) \\
\texttt{uFit\_inf}         & $n\times 1$ & $u_{\hat y_i}$ from \texttt{covBeta\_inf} \\
\texttt{uTotal\_chi2scaled}& $n\times 1$ & Per-point $u$, Option~A \\
\texttt{uTotal\_wrmsAdded} & $n\times 1$ & Per-point $u$, Option~B \\
\texttt{uBar\_chi2scaled}  & scalar      & $\bar{u}$ summary, Option~A (fit-space) \\
\texttt{uBar\_wrmsAdded}   & scalar      & $\bar{u}$ summary, Option~B (fit-space) \\
\texttt{uEff}              & $n\times 1$ & $u_{\mathrm{eff},i}$ \\
\texttt{X}, \texttt{uX}   & $n\times p$ & Design matrix and its uncertainty \\
\texttt{Xhat}              & $n\times p$ & Corrected design matrix (TLS; $=\mathbf{X}$ for WLS) \\
\texttt{x\_valid}, \texttt{y\_valid} & $n\times 1$ & Data after NaN removal \\
\texttt{method}            & \texttt{str}& \mstr{WLS} or \mstr{TLS} \\
\texttt{model\_str}        & \texttt{str}& Model string as passed \emph{(note: snake\_case)} \\
\texttt{params}            & \texttt{list[str]} & Parameter names as passed \\
\texttt{n}, \texttt{p}    & \texttt{int}& Observations used; parameters \\
\bottomrule
\end{tabular}

% ─────────────────────────────────────────────────────────────────────────────
\subsection{Examples}
% ─────────────────────────────────────────────────────────────────────────────

\subsubsection*{Example 1: Steinhart--Hart thermistor calibration}

The Steinhart--Hart model relates thermistor resistance $R$ to temperature
$T$ via
\begin{equation}
  \frac{1}{T} = A + B\ln R + D(\ln R)^3,
\end{equation}
so the fit is performed in $1/T$ space and a transform function converts the
result back to display units ($^\circ$C).  The TLS call additionally
propagates resistance uncertainty $u(R)$ through the design matrix.

\begin{lstlisting}[style=matlab, caption={Steinhart--Hart fit, MATLAB.}]
% --- Data ---
R         = [...];             % resistance [Ohm]
T         = [...] + 273.15;   % temperature [K]
uR        = [...];             % standard uncertainty in R [Ohm]
uT        = [...];             % standard uncertainty in T [K]
uOneOverT = uT ./ T.^2;       % propagated uncertainty in 1/T [1/K]

% --- Model ---
modelStr = 'A + B*log(x) + D*(log(x))^3';
params   = ["A", "B", "D"];

% WLS: uncertainty in y (1/T) only
wls = LinearRegression(R, 1./T, uOneOverT, modelStr, params, ...
    'transformFcn', @(v) 1./v - 273.15, ...
    'xLabel', 'R ($\Omega$)', 'yLabel', 'T ($^\circ$C)', ...
    'residScale', 1000, 'residUnit', 'mK', ...
    'plotTitle', 'Steinhart--Hart WLS');

% TLS: uncertainty in both R and 1/T
tls = LinearRegression(R, 1./T, uR, uOneOverT, modelStr, params, ...
    'transformFcn', @(v) 1./v - 273.15, ...
    'xLabel', 'R ($\Omega$)', 'yLabel', 'T ($^\circ$C)', ...
    'residScale', 1000, 'residUnit', 'mK', ...
    'plotTitle', 'Steinhart--Hart TLS');

% --- Access results ---
fprintf('WLS  A = %.4e +/- %.4e\n', wls.beta(1), wls.uBeta(1));
fprintf('WLS  chi2r = %.4f\n',      wls.chi2r);
fprintf('WLS  uBar (Option B) = %.2f mK\n', wls.uBar_wrmsAdded * 1000);
\end{lstlisting}

\begin{lstlisting}[style=python, caption={Steinhart--Hart fit, Python.}]
import numpy as np
from linear_regression import linear_regression, make_plot

# --- Data ---
R         = np.array([...])            # resistance [Ohm]
T         = np.array([...]) + 273.15  # temperature [K]
uR        = np.array([...])            # standard uncertainty in R [Ohm]
uT        = np.array([...])            # standard uncertainty in T [K]
uOneOverT = uT / T**2                  # propagated uncertainty in 1/T [1/K]

# --- Model ---
model_str = 'A + B*log(x) + D*(log(x))^3'   # ^ auto-converted to **
params    = ['A', 'B', 'D']

# WLS: uncertainty in y (1/T) only
wls = linear_regression(
    R, 1.0/T, uOneOverT, model_str, params,
    transform_fcn = lambda v: 1.0/v - 273.15,
    x_label    = r'R ($\Omega$)',
    y_label    = r'T ($^\circ$C)',
    resid_scale = 1000,
    resid_unit  = 'mK',
    plot_title  = 'Steinhart-Hart WLS')

# TLS: uncertainty in both R and 1/T
tls = linear_regression(
    R, 1.0/T, uR, uOneOverT, model_str, params,
    transform_fcn = lambda v: 1.0/v - 273.15,
    x_label    = r'R ($\Omega$)',
    y_label    = r'T ($^\circ$C)',
    resid_scale = 1000,
    resid_unit  = 'mK',
    plot_title  = 'Steinhart-Hart TLS')

# --- Access results ---
print(f"WLS  A = {wls.beta[0]:.4e} +/- {wls.uBeta[0]:.4e}")
print(f"WLS  chi2r = {wls.chi2r:.4f}")
print(f"WLS  uBar (Option B) = {wls.uBar_wrmsAdded * 1000:.2f} mK")

# --- Figures (separate call in Python) ---
for res, name in ((wls, 'wls'), (tls, 'tls')):
    fig = make_plot(
        res,
        transform_fcn = lambda v: 1.0/v - 273.15,
        x_label     = r'R ($\Omega$)',
        y_label     = r'T ($^\circ$C)',
        resid_scale = 1000,
        resid_unit  = 'mK',
        plot_title  = f'Steinhart-Hart {name.upper()}')
    fig.savefig(f'steinhart_hart_{name}.png', dpi=150)
\end{lstlisting}

% ── ── ── ── ── ── ── ── ── ── ── ── ── ── ── ── ── ── ── ── ──

\subsubsection*{Example 2: Quadratic polynomial}

A simple quadratic $y = A + Bx + Cx^2$ with uniform $y$-uncertainties.
No transform function is needed because the fit space and display space
are the same.

\begin{lstlisting}[style=matlab, caption={Quadratic polynomial fit, MATLAB.}]
% --- Model ---
modelStr = 'A + B*x + C*x^2';
params   = ["A", "B", "C"];

% WLS
result = LinearRegression(x, y, uY, modelStr, params, ...
    'xLabel', 'x', 'yLabel', 'y', ...
    'plotTitle', 'Quadratic fit');

% --- Access results ---
fprintf('A = %.4g,  B = %.4g,  C = %.4g\n', ...
    result.beta(1), result.beta(2), result.beta(3));
fprintf('Reduced chi-squared: %.4f\n', result.chi2r);
\end{lstlisting}

\begin{lstlisting}[style=python, caption={Quadratic polynomial fit, Python.}]
# --- Model ---
model_str = 'A + B*x + C*x**2'   # or 'A + B*x + C*x^2' -- both work
params    = ['A', 'B', 'C']

# WLS
result = linear_regression(x, y, uY, model_str, params,
    x_label    = 'x',
    y_label    = 'y',
    plot_title = 'Quadratic fit')

# --- Access results ---
A, B, C = result.beta
print(f"A = {A:.4g},  B = {B:.4g},  C = {C:.4g}")
print(f"Reduced chi-squared: {result.chi2r:.4f}")
\end{lstlisting}