\documentclass[twoside]{MATH77}
\usepackage{multicol}
\usepackage[fleqn,reqno,centertags]{amsmath}
\begin{document}
\hyphenation{SCPVAL}
\begmath 11.1 Polynomial Least-Squares Curve Fit

\silentfootnote{$^\copyright$1997 Calif. Inst. of Technology, \thisyear \ Math \`a la Carte, Inc.}

\subsection{Purpose}

This subroutine determines a univariate polynomial that fits a given discrete
set of data in the sense of minimizing the weighted sum of squares of residuals.
The fitted polynomial can be constrained to match some data points exactly
by appropriate setting of the a priori standard deviations of the data
errors.

Auxiliary subroutines, described in Chapter~11.2, may be used to evaluate,
differentiate, or integrate a polynomial produced by this fitting procedure.

\subsection{Usage}

\subsubsection{Program Prototype, Single Precision}
\begin{description}
\item[INTEGER]  \ {\bf M, NMAX, NDEG}

\item[REAL]  \ {\bf X}($\geq $M),{\bf Y}($\geq $M),{\bf SD}($\geq $M~or~=1),%
{\bf P}($\geq $NMAX+3), {\bf SIGFAC, W}($\geq$(NMAX+3)$\times $(NMAX+3))

\item[LOGICAL]  \ {\bf SEEKN, COMTRN, CHBBAS}
\end{description}
Assign values to M, X(), Y(), SD(), NMAX, SEEKN, COMTRN and CHBBAS. If
COMTRN = .FALSE. values must also be assigned to P(1) and P(2).%
\begin{center}
\fbox{\begin{tabular}{@{\bf }c}
CALL SPFIT (M, X, Y, SD, NMAX, SEEKN,\\
COMTRN, CHBBAS, P, NDEG,SIGFAC, W)\\
\end{tabular}}
\end{center}
Following the CALL to SPFIT one may wish to use SCPVAL or SMPVAL to evaluate
the fitted polynomial at specific points, SCPDRV or SMPDRV to differentiate
the polynomial, or SCPINT or SMPINT to integrate the polynomial. See Chapter
11.2 for descriptions of these subprograms.

\subsubsection{Argument Definitions}
\begin{description}
\item[M]  \ [in] Number of data points.

\item[X()]  \ [in] Array of values of the independent variable. Values need
not be ordered and may be repeated.

\item[Y()]  \ [in] Array of values of the dependent variable indexed to
correspond to the X() values.

\item[SD()]  \ [in] If SD($1)>0.$, each SD($i$) must be positive and will be
assumed to be the user's a priori estimate of the standard deviation of the
uncertainty (observational errors, etc.) in the corresponding data value
Y($i$).

If SD($1)<0.$, $|$SD(1)$|$ will be used as the a priori standard deviation
of each data value Y($i$). In this case the array SD() can be dimensioned
SD(1) since locations following SD(1) will not be used.

\item[NMAX]  \ [in] Specifies the maximum degree polynomial to be
considered. Require NMAX\ $\geq 0$. If NMAX $>$ M $-$ 1 the subroutine will
function as though NMAX = M $-$ 1.

\item[SEEKN]  \ [in] If SEEKN = .TRUE. the subroutine will determine the
optimum value of NDEG $\leq $ NMAX in the sense described below in Section D.

If SEEKN = .FALSE. the subroutine will set NDEG = NMAX unless this is a
singular or nearly singular problem, in which case NDEG will be reduced as
necessary.

\item[COMTRN]  \ [in] If COMTRN = .TRUE. the subroutine will determine the
transformation parameters P(1) and P(2) as described in Section D.

If COMTRN = .FALSE. the subroutine will use the values of P(1) and P(2) that
have been set by the user.

\item[CHBBAS]  \ [in] If CHBBAS = .TRUE. the subroutine will use the
Chebyshev basis. If CHBBAS = .FALSE. the subroutine will use the monomial
basis. The Chebyshev basis is recommended for better numerical stability,
when NMAX is greater than about six.

\item[P()]  \ [inout] Parameter vector defining a polynomial according to
Eqs.\,(1) and (2) or (3) and (4). Values of P(1) and P(2) may either be set
by the user, or by the subroutine, depending upon the value of COMTRN.
Values of P($j$+3), $j$ = 0, ..., NMAX will be computed by the subroutine. If
NDEG $<$ NMAX, the subroutine will set P($i$+3) = 0, for $i$ = NDEG + 1, ...,
NMAX.

\item[NDEG]  \ [out] Set by the subroutine to indicate the degree of the
fitted polynomial. On return NDEG will satisfy 0 $\leq $ NDEG $\leq $ NMAX
if the computation was successful and will be set to $-$1 if an error occurred.

\item[SIGFAC]  \ [out] Number computed by the subroutine. See discussion of
SIGFAC in Sections C and D following.

\item[W()]  \ [scratch] Array of working storage.
\end{description}
\subsubsection{Modifications for Double Precision}

For double precision usage change the REAL type statements to DOUBLE
PRECISION and change the subroutine name from SPFIT to DPFIT. Change the
names of the auxiliary subroutines mentioned above to DCPVAL, DMPVAL,
DCPDRV, DMPDRV, DCPINT, and DMPINT, respectively.

\subsection{Examples and Remarks}

\subsubsection{Example}

Given a set of 12~data pairs $(x_i$, $y_i)$, compute the
weighted least-squares polynomial fit to these data letting the subroutine
determine the preferred polynomial degree not to exceed~8. After computing
the least-squares polynomial $p(x)$, compute and tabulate the quantities $%
x_i $, $y_i$, $p(x_i)$, and $r_i = y_i - p(x_i).$

The program DRSPFIT carries out this computation using subroutines SPFIT and
SCPVAL. The output is shown in ODSPFIT. Note that the value of NDEG selected
by the subroutine is~7.

\subsubsection{Transformation of the Independent Variable}

For the best numerical accuracy, both in determining the best fitting
polynomial and in the later evaluation of the polynomial, it is generally
advisable to use transformation parameters P(1) and P(2) that cause the
transformed variable $u$ of Eqs.\,(1) or (3) to range from $-$1 to +1. This
condition can be obtained by setting COMTRN = .TRUE.

When fitting with polynomials of degree higher than about six, the Chebyshev
basis generally gives better numerical accuracy than the monomial basis. To
obtain the potential numerical advantages of the Chebyshev basis it is
essential to cause the transformed variable, $u$, to range from $-1$ to +1.

To force an identity transformation $(i.e.$, effectively no transformation)
of the independent variable, set COMTRN = .FALSE., P(1) = 0., and P(2) = 1.

\subsubsection{Interpretation of SIGFAC}

If the user has set all SD($i$) = 1.0, or equivalently, set SD(1) = $-$1.0,
SIGFAC can be interpreted as an a posteriori estimate of the standard
deviation of the error in each Y($i$) value. More generally, whatever values
the user has assigned to the SD($i$)'s, the a posteriori estimate of the
standard deviation of the error in Y($i$) is SIGFAC $\times $ SD($i$).

\subsubsection{Equality Constraints}

If the user wishes the fitted polynomial to agree to machine accuracy with
one or more of the data points, this can be accomplished by setting the SD()
value for such points much smaller than for the other points.

Let $\epsilon $ denote the machine precision, $i.e.$, $\epsilon =$
R1MACH(4) in single precision or $\epsilon =$ D1MACH(4) in double
precision, see Chapter~19.1.  For some value of $c$
in the range 0.5 $\leq c \leq 0.75$, we suggest setting SD($i) = \epsilon
^c$ for the points at which an exact fit is desired and SD($i) = 1.0$ for
all other points. For example if $\epsilon = 10^{-8}$ set SD($i$) in the
range $10^{-4}$ through $10^{-6}$ for the exact fit points.

Using $c < 0.5$ may not produce sufficiently small residuals at the desired
points of exact fit, whereas setting $c > 0.75$ may trip the
near-singularity test in these subroutines leading to an unwanted alteration
of the problem.

Note that it is not mathematically reasonable to attempt to force exact fits
at more than NMAX$+1$ points.

\subsection{Functional Description}

The user provides data $(x_i$, $y_i$, $s_i)$, for $i=1$, ..., M, and an
integer NMAX. If SEEKN = .FALSE. and the problem is not rank-deficient or
extremely ill-conditioned, the subroutine simply determines the polynomial $%
p_n(x)$ of degree $n=$ NMAX that minimizes the following weighted sum of
squares of residuals:
\begin{equation*}
\rho _n^2=\sum_{i=1}^M\left[ \frac{y_i-p_n(x_i)}{s_i}\right] ^2
\end{equation*}
The subroutine also computes%
\begin{equation*}
\sigma _n=\left[ \frac{\rho _n^2}{\max (1,M-n-1)}\right] ^{1/2}
\end{equation*}
and returns this value as SIGFAC.

The subroutine makes use of subroutines SROTMG and SROTM (or DROTMG and
DROTM) from the BLAS1 (See Chapter~6.3.) to implement a sequential
Modified Givens orthogonal transformation method of reducing the matrix of
the least-squares problem to triangular form \cite{Lawson:1974:SLS}, and
then solves the triangular system.  This method has excellent numerical
stability.

After triangularizing the matrix of the problem for degree NMAX, the
quantities $\rho ^2_n$ and $\sigma ^2_n$ may be computed inexpensively for
all degrees $n$ from zero through NMAX. Thus, if the user sets SEEKN =
.TRUE. the subroutine computes these values of $\rho ^2_n$ and computes
\begin{equation*}
\text{CMIN} = \min \{\sigma ^2_n : 0 \leq n \leq \text{NMAX}\}.
\end{equation*}
The quantity $\rho _n$ is a strictly nonincreasing function of $n$ whereas $%
\sigma _n$ typically decreases initially with increasing $n$ but then
oscillates and slowly increases as $n$ becomes so large that the polynomial $%
p_n(x)$ is fitting noise rather than a true polynomial signal in the data.

The subroutine sets NDEG to be the smallest value of $n$ for which $\sigma
^2_n \leq 1.01 \times $CMIN, and it sets SIGFAC $= \sigma _{NDEG}$. This
method of selecting NDEG is generally satisfactory when the ratio NMAX/M is
reasonably small, say less than $1/3$. If this ratio is larger there is a
tendency for the NDEG selected by this method to be larger than one might
pick if one viewed graphs of solution polynomials of different degrees.

Whatever the setting of SEEKN, if the problem for degree NMAX is
rank-deficient, or very ill-conditioned, the subroutine will restrict
consideration to lower degrees for which the problem is not extremely
ill-conditioned. In particular this reduction of degree will certainly occur
if NMAX$+1 > $ M.

\subparagraph{Parameterization of the polynomial, $p_n(x)$}

In this subroutine an $n^{th}$ degree polynomial $p_n(x)$ is represented in
one of two special parametric forms (monomial basis or Chebyshev basis),
both of which include the specification of a linear transformation of the
independent variable in the parameterization. The parameters defining an $%
n^{th}$ degree polynomial $p_n(x)$ are stored in the array P($i$), $i=1$, ...,
$n+3.$

Using the monomial basis these parameters define the polynomial $p_n(x)$ as
follows:
\begin{align}
\label{O1}u&=\frac{x-P(1)}{P(2)}\\
\label{O2}p_n(x)&=\sum_{i=0}^nP(i+3)u^i
\end{align}
whereas, if the Chebyshev basis is requested, $p_n(x)$ is defined as:
\begin{align}
\label{O3}u&=\frac{x-P(1)}{P(2)}\\
\label{O4}p_n(x)&=\sum_{i=0}^nP(i+3)T_i(u)
\end{align}
where $T_i(u)$ denotes the Chebyshev polynomial of degree $i$. The Chebyshev
polynomials may be defined as follows:%
\begin{gather*}
T_0(u)=1,\quad T_1(u)=u,\quad T_2(u)=2u^2-1\vspace{-4pt}\\
T_i(u)=2u\ T_{i-1}(u)-T_{i-2}(u)\quad i=3,4,...
\end{gather*}
The parameters P(1) and P(2) that occur in the transformation formula (1) or
(3) may be set by the user (COMTRN = .FALSE.) or else computed by this
subroutine (COMTRN = .TRUE.). In the latter case P(1) and P(2) will be
computed as\vspace{-2pt}%
\begin{align*}
\text{P}(1)&=(x_{\max }+x_{\min })/2\\
\text{P}(2)&=(x_{\max }-x_{\min })/2
\end{align*}
where $x_{\max }$ and $x_{\min }$ are respectively the maximum and minimum
values of $x$ in the data array X($i$), $i=1$, ..., M. This causes the
transformed variable $u$ of Eq.\,(1) or (3) to range from $-1$ to +1 as $x$
ranges from $x_{\min }$ to $x_{\max }.$

\bibliography{math77}
\bibliographystyle{math77}

\subsection{Error Procedures and Restrictions}

The use of the Chebyshev basis is effective in improving the conditioning of
the problem only if a transformation of the variable is used that maps the
interval of interest for both fitting and evaluation to an interval
approximately coincident with [$-1.$,~1.].

The automatic procedure for selecting the degree of the fitted polynomial,
used when SEEKN = .TRUE., tends to select the degree somewhat higher in some
cases than one would choose by looking at plots of fits of different degrees.

This subroutine will issue an error message using the error processor
in Chapter 19.2 at level 0, set NDEG $= -1$, and return if
any of the following erroneous conditions exists:
\begin{itemize}
\item[1.]  SD($1)=0$, or M $\leq $ 0, or NMAX $<$ 0.

\item[2.]  COMTRN = .FALSE. and P(2) = 0.

\item[3.]  SD($1)>0$. and some SD($i) \leq 0$. for $1< i \leq $ M.
\end{itemize}
\subsection{Supporting Information}

The source language is Fortran~77.

\begin{tabular}{@{\bf}l@{\hspace{5pt}}l}
\bf Entry & \hspace{.35in} {\bf Required Files}\vspace{2pt} \\
DPFIT & \parbox[t]{2.7in}{\hyphenpenalty10000 \raggedright
AMACH, DERM1, DERV1, DPFIT, DROTM, DROTMG, ERFIN, ERMSG, IERM1, IERV1\rule[-5pt]{0pt}{8pt}}\\
SPFIT & \parbox[t]{2.7in}{\hyphenpenalty10000 \raggedright
AMACH, ERFIN, ERMSG, IERM1, IERV1, SERM1, SERV1, SPFIT, SROTM, SROTMG}\\
\end{tabular}

Initially designed and programmed by C.~L.~Lawson, JPL, 1970. Modified by
Lawson and S.~Y.~Chiu, 1984, to use Modified Givens rather than Householder
transformations to reduce the storage requirement for W() and the complexity
of the program. Also, adapted the code to Fortran~77.


\begcodenp
\lstset{language=[77]Fortran,showstringspaces=false}
\lstset{xleftmargin=.8in}

\centerline{\bf \large DRSPFIT}\vspace{0pt}
\lstinputlisting{\codeloc{spfit}}
\newpage

\vspace{30pt}\centerline{\bf \large ODSPFIT}\vspace{10pt}
\lstset{language={}}
\lstinputlisting{\outputloc{spfit}}
\end{document}

.