\chapter{Fitting}

\section{Errors}

This section describes a QDP file that contains errors
and how to control the plotting of those errors.
The next two sections describe how to define a model,
find the best fitting parameter values,
and then estimate the uncertainties on the parameter values.
Although the examples will be based on the QDP file containing errors,
it is possible (and sometimes better) to fit data without errors.

You should now create a {\tt DEMO1.QDP} file that contains
the following:
\begin{verbatim}
 READ Serr 1 2
 LAbel X Time
 LAbel Y Distance
  1.0  .25   1.24  .3
  1.5  .25   1.86  .3
  2.0  .25   3.76  .3
  4.0 1.75  16.43  .3
  7.0 1.25  49.06  .3
\end{verbatim}
The first line in this file is not a PLT command but rather
a QDP command.
The QDP {\tt READ} command is used by QDP
to tell PLT which vectors contain errors.
In this case the READ command tells QDP that vectors 1 and 2
will have symmetric errors;
hence, columns 1 and 2 contain data and errors for vector 1,
and columns 3 and 4 contain data and errors for vector 2.
Following the QDP command are two PLT commands that will be passed
to the PLT program and executed before the graph is drawn.
Including PLT commands in the QDP file provides a way to override
built-in defaults and/or to add labels to the graph.
Data lines occur after all the command lines.
To read and plot this file, use
\begin{verbatim}
$ QDP DEMO1
(enter device type)
PLT>
\end{verbatim}
\begin{figure}
   \vspace{10.6cm}
   \special{psfile=fig41.ps hscale=70 vscale=40}
   \caption{The default appearance of the {\tt DEMO1.QDP} file.}
\end{figure}
The graph which looks like Figure 4.1 should now appear.
When you plot data containing errors,
the error attribute is {\tt ON},
and the errors plotted.
As described in Section 3.3, the line will no longer appear
connecting the data.
If you want to see that line,
then should use \verb@LIne ON@ to explicitly switch on the line.
To suppress plotting of the errors use,
\begin{verbatim}
PLT> Error OFf
PLT> P
\end{verbatim}
This disables plotting of the error bars. 
Once again the line connecting the points appears.
You should now enter
\begin{verbatim}
PLT> LIne Step
PLT> P
\end{verbatim}
to produce a stepped-line plot.
Using \, {\tt Error~ON} \, at this time will cause
both the stepped line and the errors to be plotted.
This is because the \, {\tt LIne~Step} \, command sets an internal flag
that a line should be plotted
and the \, {\tt Error~ON} \, command sets another flag to plot errors.
The command \, {\tt LIne~OFf} will turn off the plotting of the (stepped) line.
The plotting of a line, errors and markers can all be turned on or off
for each vector independently.
See the {\tt LIne}, {\tt Error} and {\tt MArker} commands,
in Appendix~\ref{pltcommands}, {\bf PLT Command Summary},
for more information.

\section{Fitting}

This section requires the {\tt DEMO1.QDP} file
described in the previous section.
Before you can fit data, you must first define a model.
First, read in the data and define a constant model with
\begin{verbatim}
$ QDP DEMO1
(enter device type)
PLT> MOdel CONS
\end{verbatim}
At this point, you will be prompted for the default initial value
for the constant.
Enter \fbox{Return} to use the default,
and at the \verb@ PLT> @ prompt, type {\tt Fit}.\footnote{
The {\it FIT} subroutine minimizes $\chi^2$ using a modified version
of Bevington's {\it CURFIT} subroutine.
Bevington's book
{\it Data Reduction and Error Analysis for the Physical Sciences},
published by McGraw-Hill in 1969,
is an excellent introduction to statistics.
Anyone interested in a detailed understanding of how CURFIT works
should consult this book.}
When {\tt Fit} runs, it first tells you which plot group is being fitted
and the range over which data are being fitted.
It is important to realize that if you have used \, {\tt R X} \,
to rescale the $x$-axis so that some points are outside the range plotted,
then these points would not be included in the fit.
You cannot exclude points using the \, {\tt R~Y} \, command.
(This is intended to prevent cheating.)
You will next see the message ``Fitting 5 points in a band of 5''.
This informs you that there are 5 points in the current $x$-range.
In order to execute faster,
the FIT routine resets the minimum and maximum of the array,
to achieve the smallest range possible that includes all points
in the $x$-range,
and so the ``band of 5'' output informs you
how big this minimum range is.
Next the FIT routine prints the current parameter values
(1.00000 in this case).
The program then prints the current value
of the weighted variance {\tt W-VAR}.
If you have errors on your data,
the weighted variance is $\chi^2$;
for no errors, {\tt W-VAR} is just the variance.
The number in {\tt ()} is $\log \lambda$
and is for the expert's use.

The {\it CURFIT} routine will terminate when the change in $\chi^2$
or, for an unweighted fit the relative change in the variance, is
less than 0.05.
If this condition has not been met after 10 iterations,
you will be prompted ``{\tt Continue fitting?~(n)}''.
Answer {\tt Y} to continue or {\tt N} to terminate.
If you are are fitting in background or batch mode,
then you should always leave a blank line after the \verb@Fit@ command.
Thus if \verb@Fit@ does not terminate, the ``{\tt Continue fitting?~(n)}''
question will read the blank line with the default answer of ``no'' and
terminate.  If the fit does terminate, PLT will see the blank line,
and ignore it.
When {\it CURFIT} terminates,
the current parameter values are again printed
and this model is drawn on the current plot.

For the above, the total variance is 18323 and hence a CONS does not
look like a very good model.  Let's try a more complicated model with
\begin{verbatim}
PLT> MOdel CO LI QU
\end{verbatim}
to include constant, linear and quadratic components.
Again you can default on all the initial values.
When you type {\tt Fit},
you should find that {\tt W-VAR} has decreased to 4.23.
This is clearly a better fit.

To generate a list of all possible built-in components,
use the \, {\tt MOdel~?} \, command.
To obtain a description of what a component does,
use the \, {\tt HElp~MOdel} \, command followed by the component name.
If you can not construct your model from the built-in components,
then you can create additional components.
A COD file can be used to define a sophisticated component.
COD files are ASCII text files that contain functions
written in a Forth-like computer language.
Chapter 7 and Appendix~\ref{codcommands} describe COD in some detail.
If your component is too complicated for COD,
or you don't like using COD,
then you can create a Fortran function that can be used as a new component.
The next chapter will describe how to create this function {\tt UFNY},
and the supporting routines to replace the built-in \verb@DEMO@
component.

It is possible to save the current model to a disk file
using the {\tt WModel} command.
For example,
\begin{verbatim}
PLT> WModel DEMO1
\end{verbatim}
will create a {\tt DEMO1.MOD} file.
To read this model back into PLT use the command \, {\tt MOdel~@DEMO1}.
Model files can be printed out
to make a hardcopy of the current parameter values.
If you do not enter a file
name with the {\tt WModel} command,
then the model is written to your current terminal screen.

\section{Parameter uncertainties}

The \, {\tt Uncertainty} \, command
can be used to estimate the uncertainties in the parameter values.\footnote
{Uncertainties in parameter values are estimated
using the method described in
``Parameter Estimation in X-ray Astronomy'',
by M. Lampton, B. Margon, and S. Bowyer,
in {\it The Astrophysical Journal} (1976) Vol. 208, p. 177.}
To try this, you should first fit the data in the {\tt DEMO1.QDP} file
to a \, {\tt CO~LI~QU } \, model,
as described in the previous section.
Now enter the command
\begin{verbatim}
PLT> Uncertain 1
\end{verbatim}
The program will now change the value of parameter 1 by a small amount
and recompute $\chi^2$.
At each step,
the delta parameter value and the $\Delta$$\chi^2$ are printed out.
For complicated models it may take many steps
before the desired value of $\Delta$$\chi^2$ is found.
The program considers both positive and negative delta parameter values.
The default value of $\Delta$$\chi^2$ is 2.7
which corresponds to the 90\% confidence range for a single parameter.

With the {\tt DEMO1.QDP} file you will find
that both parameter 1 and 2 are consistent with zero.
To see whether they can be eliminated, try the following:
\begin{verbatim}
PLT> Newpar 1,0,-1
PLT> Newpar 2,0,-1
PLT> Fit
\end{verbatim}
The first command resets both the {\tt VAL} and {\tt SIG} terms
of parameter 1 to be 0.0 and -1.0, respectively.
A {\tt SIG} of -1.0 means that the parameter is frozen
and hence not allowed to change.
Note: The command \, {\tt Newpar~1,,-1} \,
would have frozen parameter 1 at its current value.
The second line freezes the value of parameter 2 to be zero.
The results of the {\tt Fit} reveal that $\chi^2$
has increased by 0.55 and the $F$-statistic
or a likelihood ratio tells us
that these two components were not required by the model.

The \, {\tt Uncertainty} \, command is fairly robust
but on occasion can have difficulties.
Sometimes, {\tt Uncertainty} will find a new minimum value of $\chi^2$.
This causes the search to be stopped
and the parameter values to be reset.
At this point, you should re-issue the {\tt Fit} command
to locate precisely the new minimum.
Sometimes, {\tt Uncertainty} will be unable
to locate the requested value of $\Delta$$\chi^2$ after 10 tries.
At this point the message \, {\tt UNCERT--Give up.} \, is printed.
It will be up to you to decide whether the error has been correctly calculated.
Finally, the \, {\tt Uncertainty} \, command
uses the {\tt SIG} value to estimate the location of the error.
If this number is greatly in error,
then \verb@Uncertainty@ will be starting its search in the wrong place.
If this occurs, then it is sometimes possible to adjust \verb@SIG@
to be a more accurate estimate,
before issuing the \verb@Uncertainty@ command.
It is also possible that the \verb@SIG@ is inaccurate because the
true minimum has not been found and further fitting is needed.
