Physics 128L Nonlinear Model Fit examples

(*First I will invent some data.  Ignore this step for your real lab data!  My made-up data will be the sum of several components: a constant ... *) 
fakedata = Table[3.45 ,{i,1,30}];
(* ... and a Gaussian peak centered around 19... *) 
fakedata = Table[fakedata[[i]] + 5.67*E^-((i-18.78)/.978)^2,{i,1,30}];
(* ... and, what the heck, some Gaussian random noise. I'll make the noise worse on the right, that's the sort of thing that might happen.*)
fakedata = Table[fakedata[[i]] +RandomReal[NormalDistribution[0,0.1+i*i/500]],{i,1,30}];
myErrors =Table[0.1+i*i/500,{i,1,30}]

(* and I'll build that into a ListPlot-style table. *) 
datalist = Table[{i,fakedata[[i]]},{i,1,30}]
ListPlot[datalist,PlotRange->All]
(* I'll toss in the error-bar version too. *)
Needs["ErrorBarPlots`"]
datalistwitherrors = Table[{{i,fakedata[[i]]},ErrorBar[myErrors[[i]]]},{i,1,30}]
myListPlot = ErrorListPlot[datalistwitherrors,PlotRange->All]

(* OK, I generally expect this data to look like the following function: *)
dumbidea = Plot[a +  b*E^(-((x-c)/d)^2), {x,1,30}]
(* but that's not plottable.  You can't plot variables, and in the above equation a,b,c,d, are just variables.   You can only plot them if everything is numeric. *)
guessedLinePlot = Plot[1.0 +  3.0*E^(-((x-4.0)/5.0)^2),{x,1,30}]
(* that's obviously not QUITE the right function, because the "true" values of the parameters are not 1,2,3,4,5.  Therefore, this attempt at the function did NOT give me a function that describes where the data is:*)
Show[myListPlot,guessedLinePlot] 
(* In principle, one thing you could do is KEEP GUESSING.   Improve your guesses---for example, make the Gaussian narrower and move its center to the right---until you find agreement.   That's a perfectly fine thing to do, but tedious. *)

(* "Fit" functions are built to automate the trial-and-error process of trying-and-retrying these functions.    The simple function "Fit" or |"FindFit" is easy to use, but it only works with  *multipliers.*   It can't find parameters that go into exponents, or trig functions, or whatever.  You give it a list of terms, it figures how to add them together.*) 
myDumbFit= Fit[datalist,{1, E^(-((x-4.0)/5.0)^2)},x]
Show[myListPlot,Plot[myDumbFit,{x,1,30}]]
(* Notice that it made a pretty good guess at the constant-offset, but it had no way to manipulate the width or position of the bump.  In this case it actually gave us a negative bump.*)

(* To do anything else we need a NONLINEAR fit, with the power to manipulate any parameter anywhere.  *)
myNLM = NonlinearModelFit[datalist,a +  b*E^(-((x-c)/d)^2),{a,b,c,d},x]
myNLMPlot = Plot[myNLM[x],{x,1,30},PlotRange->All]

Show[myNLMPlot,myListPlot]
myNLM["ParameterTable"]
(* amazing!  It worked.   Sometimes it doesn't work --- a common problem is that the automated "guesswork" process has to start somewhere, and improve gradually.  (The default is to start ALL of the parameters at 1.000.)   If the first few guesses aren't going in the right direction, it gives up ... sometimes by failing, sometimes by finding a god-awful fit (this is called "getting stuck in a local minimum") Always plot your results to make sure your fit is reasonable*)

(* the NonlinearModelFit will almost always work IF you have a good guess at the parameters.   Tell the fitter to start with a reasonable ballpark guess, and it'll do all of the optimization smoothly.   In this case, I can provide some OK guesses by just eyeballing the function.*)

Plot[4 +  5*E^(-((x-20)/1)^2),{x,1,30},PlotRange->All]
myNLMGood = NonlinearModelFit[datalist,a +  b*E^(-((x-c)/d)^2),{{a,4},{b,5},{c,20},{d,1}},x]
Show[Plot[myNLMGood[x],{x,1,30},PlotRange->All],myListPlot]
myNLMGood["ParameterTable"]
(* Perfect!  Notice that the fitter has gotten "the right answers" in this case---it has found values consistent (i.e. within about an error bar) with the "truth" that was responsible for generating this data in the first place.*)

(* The next question is about error bars.   What do you do if your data points have errors on them?   How do you tell NonLinearFit about the error bars?   (Notice that I did not feed the error-bar data into the fit, just the points.)  It turns out that what you need to do is to *weight* the points.  A point with a large error bar E is downweighted by 1/El, so the fitter relies on it less heavily.   A point with a small error bar Es should be weighted heavily, as 1/Es.  Let's turn our list of errors into a list of weights.*)
(*I'll start with the errors that I invented for this problem.  In reality, you should use your actual errors*)

myWeights = 1/myErrors
myNLMWeighted = NonlinearModelFit[datalist,a +  b*E^(-((x-c)/d)^2),{{a,4},{b,5},{c,20},{d,1}},x,Weights->myWeights]

(* I'll compare the data and the old and new fits. *)
Show[Plot[myNLMWeighted[x],{x,1,30},PlotRange->All],Plot[myNLMGood[x],{x,1,30},PlotRange->All],myListPlot]
{myNLMWeighted["ParameterTable"],myNLMGood["ParameterTable"]}