Lsq2top

Purpose

Fits a polynomial to the top/(bottom) of data.

Synopsis

[b,resnorm,residual,options] = lsq2top(x,y,order,res,options)

Description

LSQ2TOP is an iterative least squares fitting algorithm. It is based on a weighted least squares approach where the weights are determined at each step. At initialization the weights are all 1, then a polynomial is fit through the data cloud using least squares. When fitting to the top of a data cloud, data points with a residual significantly below a defined limit (i.e. the points below the polynomial fit line) are given a small weighting. Therefore, on subsequent iterations these data points are weighted less in the fit, and the fit line moves to fit to the top of the data cloud.

Input x is the independent variable e.g. a Mx1 vector corresponding to a frequency or wavelength axis. Input y is the dependent variable e.g. a Mx1 vector corresponding to a measured spectrum. Input order is a scalar defining the order of polynomial to be fit e.g. y = P(x), and res is a scalar approximation of the fit residual e.g. noise level. Input options is discussed below. Note that the function can be used to fit to the top or bottom of a data cloud by changing trbflag in options.

The outputs are b, the regression coefficients [highest order term corresponds to b(1) and the intercept corresponds to b(end)], resnorm is the squared 2-norm of the residual, and residual is the fit residuals = y - P(x). The options ouput is the input options echoed back, the field initwt may have been modified.

Options

• display: [ 'off' | {'on'} ] governs level of display to command window.

• trbflag: [ 'top' | {'bottom'} ] top or bottom flag, tells algorithm to fit the polynomials, y = P(x), to the top or bottom of the data cloud.

• tsqlim: [ 0.99 ] limit that governs whether a data point is significantly outside the fit residual defined by input res.

• stopcrit: [1e-4 1e-4 1000 360] stopping criteria, iteration is continued until one of the stopping criterion is met: [(relative tolerance) (absolute tolerance) (maximum number of iterations) (maximum time [seconds])].

• initwt: [ ] empty or Mx1 vector of initial weights (0<=w<=1).

Algorithm

For order = 1 and fitting to the top of a data cloud, LSQ2TOP finds the vector $\mathbf {b} =\left[{\begin{matrix}b_{1}&b_{2}\\\end{matrix}}\right]$ that minimizes $\left(\mathbf {y} -\mathbf {x} b_{1}-\mathbf {1} b_{2}\right)^{T}\mathbf {W} \left(\mathbf {y} -\mathbf {x} b_{1}-\mathbf {1} b_{2}\right)$ where W is a diagonal weighting matrix whose elements are initially 1 and then are modified on each subsequent iteration.

The weighting is determined by first estimating the residuals for each data point j as $residual_{j}=\mathbf {y} _{j}-\mathbf {x} _{j}b_{1}-b_{2}$ and defining $t_{j}=residual_{j}/res$ where res is the input res. A corresponding t-statistic from a t-table is estimated using the following

tsqst = ttestp(1-options.tsqlim,5000,2);

where $t_{table}$ is tsqst. The elements of W are then given by $w_{j}={1}/{\left(0.5+t_{j}/t_{table}\right)}\;$ for data points with $t_{j}<t_{table}$, and is a 1 otherwise. Therefore, the weighting is smaller for points far below the fit line.

The procedure can be modified to fit to the bottom of a data cloud by changing options.trbflag.

See Also

baseline, baselinew, fastnnls, fastnnls_sel, med2top