Overview of Matlab Curve Fitting Toolbox

Junior Lab Technical Sta

MIT Department of Physics
Advanced Experimental Physics I & II
This quick-start guide contains instructions on how to use Matlab on Athena to t data sets you have obtained in Junior Lab experiments. No prior knowledge of Matlab is necessary and only minimal familiarity with Athena is assumed. 1. WHY MATLAB?
create a matlab directory in your home directory. You can do this by typing: % cd ~; mkdir matlab 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 myscript.m type: % emacs myscript.m & Notice the.m extension on the script. All Matlab scripts must end in.m in order to execute correctly.
3. STARTING MATLAB ON ATHENA
Curve tting is one of the most common analytical tasks you will perform during Junior Lab. There exist many commercially available software packages for data manipulation, analysis and presentation. Some com mon programs you may have used before include Matlab, Mathematica, Origin, LabVIEW and Excel. These products vary widely in their ease of use, their power, and their level of technical support and documen tation. In order to ensure some level of uniformity and common understanding of the underlying techniques and methods implemented in a specic package, we have de cided to standardize on the use of Matlab. This standard ization has the additional benet that Matlab is available on Athena and therefore ubiquitous at MIT. You will nd that it is also a very common tool used throughout Academia and Industry so it is a package well worth learning. Athena has recently added a Curve Fitting Toolbox (cftool) to its Matlab libraries. This toolbox oers ex tensive tting capability and is quite adequate for Ju nior Lab purposes. (For more information and com plete documentation on cftool, you can visit the Mathworks website: http://www.mathworks.com/access/ helpdesk/help/toolbox/curvefit/curvefit.shtml) This short guide is designed to get you started using Matlab and the Curve Fitting Toolbox and to suggest other places for further information. Syntax and practical issues notwithstanding, remem ber that the important thing is not which software package you use as is to develop a strong understanding of the underlying mathematics, detailed in many places, most notably the required Junior Lab text, see [1].
Using Matlab requires the following steps: 1. From the Athena prompt, attach the Matlab locker (you can also add this to you.cshrc.mine le to have it automatically attached at login): % add matlab 2. Open Matlab into a new window: % matlab & Once it has nished loading, it will present you with the Matlab prompt: >>. From this prompt you can execute an of the Matlab commands or run a Matlab script. To run a script, rst make sure it ends in.m and resides in your matlab directory and then simply type the name at the prompt (without the.m): >> myscript
USING MATLAB SCRIPTS 4. USING THE JUNIOR LAB FITTING TEMPLATE SCRIPT
One very powerful and very easy way to using Matlab is to use scripts. Scripts are simply text les that con tain a series of Matlab commands. The entire process of tting 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 t, make modications, etc. without having to retype all of the commands. If you plan to use Matlab scripts, its a good idea to
To accommodate quick and easy tting for Matlab beginners, we have created a script that you can use as a template for tting. This script is available on the Ju nior Lab homepage or you can copy it directly into youre matlab directory by typing % cp /mit/8.13/matlab/fittemplate.m ~/matlab
2 The basic procedure for using this tting script is outlined as follows: 1. Open the script in Emacs. 2. Modify the script such that it: (a) Loads data from your le into the vectors x and y (b) Assigns the appropriate errors (c) Contains the functional form that you want to t (d) Uses the options needed for your t (start ing point, upper and lower bounds, lin ear/nonlinear, etc.) (e) Plots gures with labels appropriate for your data 3. Save the le with a new name, keeping the.m extension 4. Open Matlab and run the script 5. Matlab will t your data, output the information relevant to the t, and plot (1) the tted curve on top of the original data and (2) the residuals. This process will be outlined in detail below in a t to a sample data set.

4.1. Non-Linear Least Squares Example
load gauss3; Then change the lines that assign values to the x and y vectors x=xpeak; y=ypeak; (This syntax has a dierent form than that in the tem plate because the gauss3 data set contains predened vectors called xpeak and ypeak which we assign to the x and y vectors respectively.)
(b) Assigns the appropriate errors
Cftool requires that you specify the weights, wi , for doing a t. The weights are simply given by

2 wi = 1/i.

In the case of gauss3, we have a constant variance of 6.25 so we can assign the weights with the commands: sig=ones(size(x))*sqrt(6.25); wgts=1./sig.^2; which creates a weights vector of the same size as the vector x with a constant value of one over 6.25. Note: Since there are several dierent ways you could assign values to the error vector sig, the script includes several dierent assignment statements that cover a few cases. Once you have chosen one and mod ied it for the t you are doing, you will need to com ment out the other lines that assign values to sigma. To comment something out, simply put a % in front of it. Latex regards any thing following a % as a com ment and does not interpret it. Youll notice that the template is heavily commented to explain what the commands are doing.
In developing this guide, we decided to use the Statistics Reference Dataset Gauss3 available from the National Institute of Standards in Technology at http://www.itl.nist.gov/div898/strd/. This is an excellent source of information regarding the testing of data analysis algorithms. The Gauss3 dataset was installed on Athena as part of the MATLAB Curve Fitting Package and lies in the default path. It consists of two poorly resolved Gaussian peaks on a decaying exponential background and must be t using using a general (nonlinear) custom model. It furthermore has normally distributed zero-mean noise with a variance of 6.25. Here we step through the process of modifying the template script to t this data. If you have not already done so, create a matlab directory and copy the template script into this directory (see above). 1. Open the script in Emacs % emacs ~/matlab/fittemplate.m & 2. Modify the script such that it: (a) Loads data from your le into the vectors x and y.
Currently the script is set up to load a dummy le called mydata.txt. Change the load command such that it will load the gauss3 data set (no.txt in this case):
(c) Contains the functional form that you want to t
The script currently contains several functional forms which you might nd useful (Linear, Gaussian, Exponential, and Lorentzian). Since we are trying to t the sum of two Gaussian s and a decaying exponen tial well need to create a custom equation for the t. We can set up the t model with the following syntax: model = fittype(a*exp(-b*x)+a1*exp(-((x b1)/c1)^2)+a2*exp(-((x-b2)/c2)^2)); Then be sure to comment out the other tting mod els so that only one is read by Latex. In general, if you need to t 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 x as the independent variable.

(d) Uses the options needed for your t (start ing point, upper and lower bounds, lin ear/nonlinear, etc.)
First you will need to decide if the t you are doing is linear or nonlinear. (Keep in mind that just because you are tting a nonlinear equation doesnt mean that you are doing a nonlinear t. The t is only nonlinear if the coecients that you are tting for appear nonlinearly in the equation). If your t is nonlinear, comment out the line that sets Value to LinearLeastSquares. If your t
is linear, comment out the line that sets Value to NonlinearLeastSquares. If you are performing a nonlinear t, your t 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 look ing at a graph of your data. The variable 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:
opts.StartPoint=[20 20]; 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.
The results of the t are stored in fresult, gof, and output: fresult contains the coecients for the t with the condence bounds. gof contains quantities describing the goodness of the t. output contains other information such as the residuals and the Jacobian
To access any of the quantities stored within these three variables, use notation of the form variable.quantity. For example, you can get the Jacobian which is stored in output by typing >> output.Jacobian If there are other quantities related to the t 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 Matlabs memory, use the command whos.
Example plot of Junior Lab data and fit results
140 T = 111.6 0.4 s 1 T2 = 147.8 0.4 s
(e) Plots gures with labels appropriate for your data
Running the script will automatically generate plots for you. Youll need to make sure these graphs display the data you want with the correct labels. Currently the script will generate one gure that is split into two sections. The top section will contain the tted curve plotted on top of your data points with errorbars. The bottom section will contain the residuals of your t. The commands title, xlabel, and ylabel label the title, x axis and y-axis of the graph respec tively. Right now the script just gives them generic titles. Make sure to modify these lines to give your graph meaningful labels with units.

80 Counts

3. Save the le with a new name, keeping the.m extension.
Be sure to save it with a dierent lename so that
the template script will be left intact for your next tting adventure. 4. 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.m). You may need to change directories if the le is stored elsewhere. Matlab has the same directory navigation commands as Athena such as ls, cd, pwd, etc. 5. Matlab will t your data, output the information relevant to the t, and plot (1) the tted curve on top of the original data and (2) the residuals. The script is currently setup to output some use ful quantities that characterize your t. These vari ables 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.

1 = 199/241 = 0.83

Probability 95%
Model: y(x)=aebx + a e((xb1) /c1)+ a e((xb2) /c2

Time (s) 250

FIG. 1: This is an example of a basic gure 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. You should check with your individual section instructor for explicit instructions on how to prepare gures for your section. The best-t parameters (at the 68.2% or 1- condence level) are: a=98.940.53, a1=100.70.8, a2=73.711.21, b=0.010950.00012, b1=111.60.4, b2=147.80.4, c1=23.30.37, c2=19.670.38.

UTILIZING MATLAB GRAPHS

The tting template script is already congured to gen
erate graphs of your data, the tted curve and the resid
uals. However, as previously mentioned you will need to
4 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. You can do this graphically with the tools provided in the graph window or by using the text command within the script. To use the text command you simply spec ify 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: text(18,5,y(x) = ae^{-bx}+a_1e^{-((x-b_1)/ c_1)^2}+a_2e^{-((x-b_2)/c_2)^2}) 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 want them (a completed graph might look like Figure 1) , 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 Print. from the File menu and make sure the correct printer is specied. If you wish to output your graph as a postscript le for use in a written summary or oral presentation, select Export. from the File menu and enter the le name you wish to use. You have a choice between many dierent le 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.

6. GETTING STARTED WITH MATLAB AT MIT
Throughout the year, and especially at the beginning of the academic year, Athena Minicourses on Matlab are oered, see web.mit.edu/afs/athena.mit.edu/astaff/ project/minidev/www/sched-year.html. If you cant make one of these (or in addition), a very nice introduction to using Matlab on Athena is available at 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.
[1] Bevington, P.R., and D.K. Robinson, Data Reduc tion and Error Analysis for the Physical Sciences, 3rd Ed.,WCB/McGraw-Hill, Boston, 2003 [2] Levenberg, K., A Method for the Solution of Certain Problems in Least Squares, Quart. Appl. Math, Vol. 2, pp. 164-168, 1944 [3] Marquardt, D., An algorithm for Least Squares Esti mation of Nonlinear Parameters, SIAM J. Appl. Math, Vol. 11, pp. 431-441, 1963 [4] Branch, M.A., T.F. Coleman, and Y. Li, A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems, SIAM Journal on Scientic Computing, Vol. 21, Number 1, pp. 1-23, 1999
to learn more about other t options that you can select. 2. Strictly speaking, the condence bounds for tted coecients are given by: C = b t S where b are the coecients of the t, t is the inverse of the Stu dents t cumulative distribution function (see Bev ington, Appendix C.6 for more details on the origin and signicance of the Student t distribution or try help tinv from inside Matlab), and S is a vec tor of the diagonal elements of the covariance ma trix of the coecient estimates, (X T X)1 s2. Here s2 = 2. Refer to Linear Least Squares from 1 Matlabs online documentation or Reference [1] for more information about X and X T. 3. The Matlab toolbox returns a goodness-of-t statis tic called SSE (Sum of Squares due to Error). This is simply the 2 statistic used frequently by physi cists and in Bevington. The other necessary statis tics is DFE (Degrees of Freedom) and is equal to the number of data points minus the number of t ted coecients. In Junior Lab, we are most often interested in the reduced chi-square (2 ) which 1 is simply given by the 2 divided by the DFE.

APPENDIX: TECHNICAL NOTES
1. Cftool has several tting algorithms at its dis posal. The most powerful tting technique is the Trust-Region discussed in [4] which is the default non-linear method. It permits the use of lower and upper bounds around the initial guesses for the tting parameters. The second method Levenberg-Marquardt is the one discussed in [1] and originally presented in [2, 3]. This method, is only slightly less powerful in that it cannot ac cept bounds on the t parameters (though it does require good initial guesses). Try typing: >> help fitoptions

Data Plotting and Curve Fitting in MATLAB

Curve Fitting

Get the file pwl.dat from the class web page. This is an ASCII text file containing two columns of numbers representing the x and y coordinates of a dataset. From MATLAB, type load pwl.dat to load the file into a matrix named pwl. Type pwl to display the matrix in text form. In order to view the data graphically, type plot(pwl(:,1), pwl(:,2), o) This plots the second column vs. the first column using circles for each data point. This is the appropriate way to display measured data from an experiment. You should not connect the points with lines because no measurements were made between the points. Incorrect way to plot measured data: plot(pwl(:,1) ,pwl(:,2)) Now that weve got the data plotted as individual points, lets label our axes: xlabel(Input Voltage [V]) ylabel(Output Current [mA]) title(Curve Fitting Exercise) It is very important to always label your axes with units! This is engineering, not mathematics we deal with physical quantities. Now lets fit this (obviously piece-wise linear) data with two straight lines. MATLAB has a handy function called polyfit that fits an n-th order polynomial to a set of data points. A line is just a 1st-order polynomial. Type help polyfit to learn more about this function. First, we need to isolate the first linear section of the data. Lets plot a subsection of the data and try to get all of the first region and none of the second region: plot(pwl(1:60,1), pwl(1:60,2), o) Oops, it looks like we plotted a bit too much of the curve. Try this again until you isolate the first region. Now we give this data to polyfit: fit1 = polyfit(pwl(1:50,1), pwl(1:50,2), 1) Look at the variable fit1. As you learned from the polyfit help file, fit1 is a vector representing the coefficients of the fitting polynomial. We can plot the line represented by these parameters by entering: plot(pwl(:,1), fit1(1)*pwl(:,1)+fit1(2)) Notice that we leave out the o parameter so the our fit is plotted as a continuous line. This is appropriate for a theoretical fit. If we want to plot the data and the fit and the same time, we can enter: plot(pwl(:,1), pwl(:,2), o, pwl(:,1), fit1(1)*pwl(:,1)+fit1(2)) If youre getting tired of typing pwl(:,1), etc., create new variables with shorter names: x = pwl(:,1); y = pwl(:,2);
Now find fit parameters for the second half of the dataset, and plot the data with both fits together: plot(x, y, o, x, fit1(1)*x+fit1(2), x, fit2(1)*x+fit2(2)) You might want to adjust the axes to zoom in on the data: axis([-20 50]) Finally, youll need to relabel the axes of the graph and provide a title.
Plotting Exponential Curves
Now load the second sample dataset: edata.dat. Like pwl, edata is a matrix of x and y coordinates. Plot the data: plot(edata(:,1), edata(:,2), o) The data seems to follow an exponential curve. Lets plot x vs. log(y): plot(edata(:,1), log10(edata(:,2)), o) [In MATLAB, log(x) = ln(x) and log10(x) = log10(x).] Its pretty clear that this data is exponential. This is a pretty poor way to represent the y-axis, however, because it only shows the exponent of the data. What would the units be? Heres the way you should plot exponential data: semilogy(edata(:,1), edata(:,2), o) The command semilogy creates a semi-logarithmic plot with nicely-labeled axes (see also semilogx and loglog). Now we can label our axes, and the y-axis units make sense: xlabel(V_{GS} [V]) ylabel(I_{D} [A]) title(Subthreshold MOSFET Behavior)
Fitting Exponential Curves
MATLAB has no command for fitting an exponential function to data. However, we can take the logarithm of the y data and fit the data with a straight line. First, lets relabel our data: x = edata(:,1); y = edata(:,2); Now, we can generate a straight-line fit to ln(y): fit = polyfit(x, log(y), 1) Now in order to reconstruct an exponential function, we have to exponentiate the fitting line y = ax+b, where fit(1) = a and fit(2) = b:

e ax +b = e b e ax

So we can plot the data with the exponential fit as: semilogy(x, y, o, x, exp(fit(2)).*exp(fit(1)*x))
The.* operator means item-by-item multiplication of two vectors. (The * operator can only be used when multiplying a scalar by a scalar or a vector by a scalar.) The exponential slope is given by fit(1). Similar techniques can be used to fit square-law data.
Other Useful MATLAB Functions
Use the help function to read about these potentially useful commands: grid hold axis gtext figure plot semilogy semilogx loglog polyfit subplot sqrt diff