\maketitle
\section{Introduction}
Curve fitting is one of the most common analytical tasks you will
perform during Junior Lab.  Students are welcome to utilize any set of
routines for curve fitting as long as the standards for reporting
results, identified in this and other Junior Lab
documents, are met.

This short guide is designed to get you started using Matlab, a
commercial product available on both Athena and the Windows PCs in
Junior Lab.  Academic versions of Matlab are also relatively
inexpensive for those students wishing to run it on their own
computers.  While Matlab is the default Junior Lab solution, some
students prefer to use alternative mathematics packages (e.g. Maple,
Mathematica or LabVIEW).  The only requirement is a solid
understanding of the underlying algorithm and mathematics, detailed in many
places, most notably the required Junior Lab text, see \cite{bevington}.  

\section{Starting Matlab}
Using Matlab on Athena\footnote{These specific instructions are for
  the Athena installation of Matlab.  Very similar (if not identical)
  operations will work on Windows platforms as well.} requires the following steps:
\begin{enumerate}
\item From the Athena prompt, attach the Matlab locker (you can also add this to you .cshrc.mine file to have it automatically attached at login):
\begin{verbatim}  % add matlab
\end{verbatim}
\item Open Matlab into a new window:
\begin{verbatim}  % matlab & 
\end{verbatim}
\end{enumerate}
Once it has finished loading, it will present you with the Matlab
prompt: {\tt >>}.  From this prompt you can execute an of the Matlab
commands or run a Matlab script.  To run a script, first make sure it
ends in {\tt .m} and resides in your {\tt matlab} directory and
then simply type the name at the prompt (without the {\tt .m}):
\begin{verbatim}  >> myscript 
\end{verbatim}

\section{Using Matlab Scripts}
One very powerful and very easy way to using Matlab is to use scripts.
Scripts are simply text files that contain a series of Matlab
commands.  The entire process of fitting will require at least a
handful of commands so it is useful to have them all in a single
script.  Once you have this script you can return to it later, repeat
your fit, make modifications, etc. without having to retype all of the
commands.

If you plan to use Matlab scripts, it's a good idea to create a
`matlab' directory in your home directory.  You can do this by typing:

\begin{verbatim}
  % cd ~; mkdir matlab
\end{verbatim}
You should save any Matlab scripts that you write in this directory.
You can use emacs (or any other text editor) to create and edit Matlab
scripts.  For example, to create a new script called {\tt myscript.m}
type:
\begin{verbatim}  % emacs myscript.m &
\end{verbatim}
Notice the ``{\tt .m}'' extension on the script.  All Matlab scripts
must end in ``{\tt .m}'' in order to execute correctly.


\section{Using the Junior Lab Fitting Template Script}
To accommodate quick and easy fitting for Matlab beginners, we have
created a script that you can use as a template for fitting.  This
script is available on the Junior Lab homepage or you can copy it (and
a sample data set upon which to perform a sample fit)
directly into you're {\tt matlab} directory by typing
\begin{verbatim}  
% cp /mit/8.13/matlab/fittemplate.m ~/matlab 
% cp /mit/8.13/matlab/gauss3.dat ~/matlab
\end{verbatim}
Before learning how to manipulate the script for your own use, try
running it by typing (within Matlab):
\begin{verbatim}
>> fittemplate
\end{verbatim}
As you can see, a Matlab scripts are is a very powerful tool with
which you will want to become very familiar!

The basic procedure for
using this fitting script is outlined as follows:
\begin{enumerate}
\item 
{\bf Open} the script in Emacs.
\item
{\bf Modify} the script such that it:
\begin{enumerate}
\item
Loads data from your file into the vectors {\tt x} and {\tt y}
\item
Assigns the appropriate errors
\item
Contains the functional form that you want to fit
\item
Uses the options needed for your fit (starting point, upper and lower
bounds, linear/nonlinear, etc.)
\item
Plots figures with labels appropriate for your data
\end{enumerate}
\item
{\bf Save} the file with a new name, keeping the .m extension
\item
Open {\bf Matlab} and run the script
\item
Matlab will {\bf fit} your data, output the information relevant to
the fit, and {\bf plot} (1) the fitted curve on top of the original
data and (2) the residuals.
\end{enumerate}
This process will be outlined in detail below in a fit to a sample data set.

\subsection{Non-Linear Least Squares Example}
In developing this guide, we have used the Statistics Reference
Dataset ``Gauss3'' available from the National Institute of Standards
in Technology at
\url{http://www.itl.nist.gov/div898/strd/nls/data/gauss3.shtml}.  Since NIST
has kindly provided certified values for this dataset, it is an
excellent test case for checking the acceptability of alternative data
fitting products.

The sample dataset, entitled ``gauss3.dat'' and residing in the
\url{http://web.mit.edu/8.13/matlab} directory, consists of two poorly
resolved Gaussian peaks on a decaying exponential background and must
be fit using using a general (nonlinear) custom model.  It has
normally distributed zero-mean noise with a variance of 6.25.

Here we step through the process of modifying the template script to
fit this data.  If you have not already done so, create a {\tt matlab}
directory and copy the template script into this directory (see
above).  Note also that the any and all parts of the ``template
script'' may be entered directly from the matlab command line.  This
is useful when diagnosing script generated errors.

\begin{enumerate}
\item {\bf Open the script in Emacs}
\begin{verbatim} % emacs ~/matlab/fittemplate.m &
\end{verbatim}

\item {\bf Modify the script such that it:}
\begin{enumerate}
\item
Loads data from your file into the vectors {\tt x} and {\tt y}.

\parindent .2 in {\footnotesize Currently the script is set up to load
  the simulated data file called `gausse3.dat' described earlier.  It
  can be changed to load whatever dataset name you require:
\begin{verbatim}
   load gauss3;
\end{verbatim}
The result of a {\tt load} command on a tab-delimited text file is a
single ``matrix'' variable with the name gauss3 (250 rows x 2 col).
The next lines in the script assign values to the {\tt x} and {\tt y}
vectors from this parent matrix.
\begin{verbatim}
   x=gauss3(:,1);
   y=gauss3(:,2);
\end{verbatim}
(Note that this syntax utilizes the entire data set in the vector
definitions.  By placing indices on either side of the colon, a subset
of the entire data file may be selected for fitting.)}
\item
Assigns the appropriate errors

{\footnotesize
Cftool requires that you specify the weights, $w_i$, for doing a fit.  The weights are simply given by 
\begin{equation}
w_i=1/\sigma_i^2.
\end{equation}
In the case of gauss3, we have a constant variance of 6.25 so we can
assign the weights with the commands:
\begin{verbatim} 
   sig=ones(size(x))*sqrt(6.25);
   wgts=1./sig.^2;
\end{verbatim}
which creates a weights vector of the same size as the vector {\tt x} with a constant value of one over 6.25. 

{\bf Note:} Since there are several different
ways you could assign values to the error vector {\tt sig}, the script
includes several different assignment statements that cover a few
cases.  Once you have chosen one and modified it for the fit you are
doing, you will need to \emph{comment out} the other lines that assign
values to sigma.  To comment something out, simply put a `{\tt \%}' in
front of it.  Latex regards any thing following a `{\tt \%}' as a
comment and does not interpret it.  You'll notice that the template is heavily commented to explain what the commands are doing. }
\item
Contains the functional form that you want to fit

\parindent .2 in{\footnotesize
The script currently contains several functional forms which you might
find useful (Linear, Gaussian, Exponential, and Lorentzian).  Since we
are trying to fit the sum of two Gaussian s and a decaying exponential
we'll need to create a custom equation for the fit.  We can set up the
fit model with the following syntax:
\begin{verbatim}
   model = fittype('a*exp(-b*x)+a1*exp(-((x-
           b1)/c1)^2)+a2*exp(-((x-b2)/c2)^2)');
\end{verbatim}
Then be sure to comment out the other fitting models so that only one
is read by Latex.  In general, if you need to fit a function that does
not appear in the script, comment out all the functions that do appear
and create a new one with the function that you need.  Be sure to use
{\tt x} as the independent variable.}
\item
Uses the options needed for your fit (starting point, upper and lower
bounds, linear/nonlinear, etc.)
 
\parindent .2 in{\footnotesize First you will need
to decide if the fit you are doing is linear or nonlinear.  (Keep in
mind that just because you are fitting a nonlinear equation doesn't
mean that you are doing a nonlinear fit.  The fit is only nonlinear if
the coefficients that you are fitting for appear nonlinearly in the
equation).  If your fit is nonlinear, comment out the line that sets
{\tt `Value'} to {\tt `LinearLeastSquares'}. If your fit is linear, comment
out the line that sets {\tt `Value'} to {\tt `NonlinearLeastSquares'}.  

If you are performing a nonlinear fit, your fit is much more likely to
be successful if you specify a reasonable starting point and
reasonable bounds.  You should be able to determine such parameters by
looking at a graph of your data.  \footnote{Cftool has several fitting algorithms at its disposal. The most
powerful fitting technique is the ``Trust-Region'' discussed in
\cite{branch} which is the default non-linear method.  It permits the
use of lower and upper bounds around the initial guesses for the
fitting parameters.  The second method ``Levenberg-Marquardt'' is the
one discussed in \cite{bevington} and originally presented in
\cite{levenberg, marquardt}.  This method, is only slightly less
powerful in that it cannot accept bounds on the fit parameters (though
it does require good initial guesses).  Try typing: {\bf help
fitoptions} to learn more about other fit options that you can select.}

The variable {\tt `opts.StartPoint'} will contain a vector with the starting
values for your parameters.  They should be assigned in alphabetical
order.  For the gauss3 data set, we have the following parameters:
\begin{verbatim}   opts.StartPoint=[100 100 40 0 100 150 20 20];\end{verbatim}

The variables opts.Lower and opts.Upper will contain vectors with the
lower and upper bounds of the parameters respectively.  These also
need to be assigned in alphabetical order.  If you decide to use these
variables, please be sure to uncomment them by removing the \% that is
currently in front of them.  }

\item
Plots figures with labels appropriate for your data

\parindent .2 in
{\footnotesize Running the script will automatically generate plots
for you.  You'll need to make sure these graphs display the data you
want with the correct labels.

Currently the script will generate one figure that is split into two
sections.  The top section will contain the fitted curve plotted on
top of your data points with errorbars.  The bottom section will
contain the residuals of your fit.

The commands {\tt `title', `xlabel'}, and {\tt `ylabel'} label the
title, x axis and y-axis of the graph respectively.  Right now the
script just gives them generic titles.  Make sure to modify these
lines to give your graph meaningful labels with units.}
\end{enumerate}
\item
{\bf Save the file with a new name, keeping the .m extension.}  \\ Be
sure to save it with a \emph{different} filename so that the template
script will be left intact for your next fitting adventure. 
\item
\parindent .2 in {\bf Open Matlab and run the script.}  \\ If the
script is saved in your matlab directory, you should be able to run it
simply by typing in the name at the Matlab prompt (without the the
{\tt .m}).  You may need to change directories if the file is stored
elsewhere.  Matlab has the same directory navigation commands as
Athena such as {\tt `ls', `cd'}, {\tt `pwd'}, etc.
\item
{\bf Matlab will fit your data, output the information relevant to the
fit, and plot (1) the fitted curve on top of the original data and (2)
the residuals.} 

The script is currently setup to output some useful
quantities that characterize your fit.  These variables will be
displayed in the Matlab window.  They script causes these particular
values to be output because they are declared with no semi-colon.  In
general, a semi-colon at the end of a Matlab command suppresses the
output of the command.

The results of the fit are stored in {\tt `fresult', `gof'}, and {\tt `output'}:\begin{itemize}
\item\underline{\tt fresult} contains the coefficients for the fit with the
confidence bounds.\footnote{Strictly speaking, the confidence bounds for fitted coefficients
are given by: $C = b \pm t \sqrt{S}$ where $b$ are the coefficients of
the fit, $t$ is the inverse of the Student's t cumulative distribution
function (see Bevington, Appendix C.6 for more details on the
origin and significance of the Student t distribution or try {\tt `help
tinv'} from inside Matlab), and $S$ is a vector of the diagonal
elements of the covariance matrix of the coefficient estimates,
$(X^TX)^{-1}s^2$.  Here $s^2 =\chi^2_{\nu-1}$.  Refer to ``Linear
Least Squares'' from Matlab's online documentation or
Reference~\cite{bevington} for more information about $X$ and $X^T$.}
\item\underline{\tt gof} contains quantities describing the goodness
of the fit.\footnote{The Matlab toolbox returns a goodness-of-fit statistic called
SSE (Sum of Squares due to Error).  This is simply the $\chi^2$
statistic used frequently by physicists and in Bevington.  The other
necessary statistics is DFE (Degrees of Freedom) and is equal to the
number of data points minus the number of fitted coefficients.  In
Junior Lab, we are most often interested in the reduced chi-square
$(\chi^2_{\nu -1})$ which is simply given by the $\chi^2$ divided by
the DFE.}
\item\underline{\tt output} contains other information such as the
residuals and the Jacobian
\end{itemize}
To access any of the quantities stored within these three variables,
use notation of the form {\tt `variable.quantity'}.  For example, you can
get the Jacobian which is stored in {\tt output} by typing
\begin{verbatim}
>> output.Jacobian
\end{verbatim}

If there are other quantities related to the fit that you need to
access, you can do so by simply typing in the name of the variable at
the Matlab command line.  For a list of the variables currently in
Matlab's memory, use the command {\tt `whos'}.

\end{enumerate}
   
\section{Generating Presentation Graphics}
 Figure \ref{fig:typical-fit-plot} demonstrates a typical graphic you
 might create for an oral presentation or written summary.
\begin{figure}[h]
  \includegraphics[width=8cm]{typical-fit-plot}
\caption{This is an example of a basic figure for Junior Lab
notebooks, presentations and written summaries. The dataset consists
of two poorly resolved gaussian peaks on an exponential background.
The noise is gaussian distributed with zero-mean and a variance of
6.25. {\bf You should check with your individual section instructor
for explicit instructions on how to prepare figures for presentation
within your section.}
%The best-fit parameters (at the 68.2\% or 1-$\sigma$ confidence level)
%are: a=98.94$\pm$0.53, a1=100.7$\pm$0.8, a2=73.71$\pm$1.21,
%b=0.01095$\pm$0.00012, b1=111.6$\pm$0.4, b2=147.8$\pm$0.4,
%c1=23.3$\pm$0.37, c2=19.67$\pm$0.38.
\label{fig:typical-fit-plot}}
\end{figure}
The fitting template script is already configured to generate graphs
of your data, the fitted curve and the residuals.  However, as
previously mentioned you will need to customize the titles and axis
labels so they are relavent and meaningful to your graph.  In
addition, you may also wish to annotate your graph with additional
information.  {\bf The fittemplate.m script has several examples of how to
automatically place properly formatted fit results within a plot.}  You
can do this graphically with the tools provided in the graph window or
by using the {\tt `text'} command within the script.  To use the {\tt
`text'} command you simply specify $x$ and $y$ coordinates and the
string you wish to appear.  The coordinates should be in the same
units of the graph to which they refer.  For example:
\begin{verbatim}
   text(18,5,'y(x) = ae^{-bx}+a_1e^{-((x-b_1)/
        c_1)^2}+a_2e^{-((x-b_2)/c_2)^2})
\end{verbatim}
would provide an appropriate label for the gauss3 graph in the lower
left-hand corner.  Notice also that Matlab can interpret latex
formating to display Greek characters, superscripts, and subscripts.

When you have your graphs just the way you (and your section
instructor!) want them (a completed graph might look like Figure
\ref{fig:typical-fit-plot}) , you can output them for use in your lab
notebooks, written summaries, and oral presentations.  If you simply
want printouts of your graphs for your notebooks, select {\tt
`Print...'} from the {\tt `File'} menu and make sure the correct
printer is specified.

If you  wish to output your  graph as a  postscript file for use  in a
written summary  or oral  presentation, select {\tt  `Export...'} from
the {\tt  `File'} menu and enter the  file name you wish  to use.  You
have a  choice between many different  file types for  the export, but
EPS (or EPS color if appropriate) will be  the most convenient choice
if you plan to use the graph in a latex document.

\section{Getting started with Matlab at MIT}

Throughout the year, and especially at the beginning of the academic
year, Athena Minicourses on Matlab are offered, see
\url{http://web.mit.edu/minidev/www/}.  If you can't make one of these
(or in addition), a very nice introduction to using Matlab on Athena
is available at \url{web.mit.edu/olh/Matlab/Matlab.html} and has links
to more extensive resources at MIT and elsewhere.  Finally, talk to
friends and classmates; many of them have a great deal of Matlab
experience.

\begin{thebibliography}{99}

\bibitem{bevington} Bevington, P.R., and D.K. Robinson, \textsl{Data
Reduction and Error Analysis for the Physical Sciences, 3rd
Ed.},WCB/McGraw-Hill, Boston, 2003
\bibitem{levenberg} Levenberg, K., ``A Method for the Solution of Certain Problems in Least Squares'', \textsl{Quart. Appl. Math}, Vol. 2, pp. 164-168, 1944
\bibitem{marquardt} Marquardt, D., ``An algorithm for Least Squares Estimation of Nonlinear Parameters'', \textsl{SIAM J. Appl. Math}, Vol. 11, pp. 431-441, 1963
\bibitem{branch} Branch, M.A., T.F. Coleman, and Y. Li, ``A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems'', \textsl{SIAM Journal on Scientific Computing}, Vol. 21, Number 1, pp. 1-23, 1999
\end{thebibliography}

\clearpage
\appendix

\section{Other Matlab Uses in Junior Lab}

\begin{figure}[h]
  \includegraphics[width=8cm]{simsingleslit}
\caption{Whenever possible, try to simulate the expected data by
plotting some functional form which incorporates the essential physics
of the problem.  Here is a single slit diffraction experiment simulation.
\label{fig:singleslit}}
\end{figure}

\begin{figure}[h]
  \includegraphics[width=8cm]{simdoubleslit}
\caption{And now, a double slit simulation...
\label{fig:doubleslit}}
\end{figure}

\begin{figure}[h]
  \includegraphics[width=8cm]{planck}
\caption{Planck Radiation Curves
\label{fig:planck}}
\end{figure}

\begin{figure}[h]
  \includegraphics[width=8cm]{mbdist}
\caption{A Matlab generated Maxwell-Boltzmann distribution
\label{fig:mbdist}}
\end{figure}


\begin{figure}[h]
  \includegraphics[width=8cm]{simexpdecay}
\caption{A Matlag generated dataset simulating an exponential decay
where the uncertainties are described by Poisson statistics.
\label{fig:expdecay}}
\end{figure}
