Mathematica workshop - worksheet 2

Organising code

Please don't send people copies of this worksheet. Instead, direct them to
http://j-star.org/mathematica_course.html
 where the latest version will always be available.
 © Jony Hudson, Imperial College London, 2006-2011

In the last worksheet we looked at some elements of Mathematica programming that come up a lot, and are a bit different from what you might be used to.

In this worksheet we'll look a little more systematically at one of the key aspects of programming in Mathematica - using functions to abstract operations. We'll do this in the context of basic data analysis tasks - the sort of thing that come up in the lab all the time.  We'll see how the Mathematica style (more generally, the functional programming style) is to build complicated systems from the bottom-up using simple ingredients. [People often argue that this bottom-up design is one of the key features of working with a functional language, as opposed to the top-down design more common with imperative languages. I think they just like to argue about stuff, but hey ...]

Plotting and fitting

Let's get warmed up by doing some simple least squares fitting. I'll start with an example and then move onto some exercises.

Here are some data I took while characterising the LIF molecule detector on the electron edm experiment. The data are a record of a photomultiplier's (PMT) output in current- and counting-modes for various levels of illumination. The data are pairs of {amplifier output (mV), count rate (kHz)}.

data = {{15,123},{17,130},{31,215},{36,253},{47,297},{58,388},{65,412},{82,510},{160,910},{250,1320},{500,2260},{870,3150},{1120,3650}};

Let's take a look at the data.

dataPlot=ListPlot[data,AxesLabel->{"amp output (mV)","count rate (kHz)"},PlotStyle->PointSize[0.02]]

We can clearly see the counting mode response becoming non-linear due to pulse pile-up in the counting electronics (don't worry if you don't really know anything about PMTs, or indeed experimental physics - it doesn't matter much for the examples). 

The Mathematica command FindFit can be used to fit a model to data. Let's model the data as a saturating exponential with an offset. I've given the fitter some reasonable guesses for the parameters.

model=a+b(1- (E^(-c x)));

fit=FindFit[data,model,{{a,100},{b,5000},{c,(10^(-3))}},x]

The best fit is returned as a set of rules.We can use these rules to transform the general model above to the best-fit model (using /.- look it up if you haven't already done so).

model/.fit

We see the replacement of the parameters with their fitted values. Let’s plot the data and fit on the same graph.

Show[
dataPlot,
Plot[model/.fit,{x,0,1200},PlotStyle->Red]
]

Looks ok. I took the data to find the calibration between count rate and voltage output in the PMT's linear region. How do we extract it? It's just the derivative of the fitted model at zero count rate.

D[model/.fit,x]/.x->0

Each mV of output is about 6000 photons/second. Note that we've got no idea what the uncertainty in this number is. Let's ignore this for the moment - we'll come back to it later.

Exercises

The code above is sufficient to fit a model to some data. In this section we'll start to abstract out the key operations. This will provide a foundation on which to perform more sophisticated data analysis tasks.

Plotting the fitted data

Write a function that takes a list of {x,y} pairs and a fitted model and plots the data and fit on the same graph. The fit plot should run over the range of x in the data list. Test your function by evaluating it with the data, model and fit from the above example (you'll need to run the cells above first).
Hint: this is quite a complicated function so you might find a Module is useful.
Hint 2: don't worry about axes labels and other bells and whistles for the moment. Unless you want to, of course!
Hint 3: you might find the functions Min and Max useful.
Hint 4: re-read the question! Lots of people misundertand it and answer a different question.

Solution

plotFit[data_,model_]:=Module[{xData,xMin,xMax},
xData=First[#1]&/@data;
xMin=Min[xData];
xMax=Max[xData];
Show[
ListPlot[data,PlotStyle->PointSize[0.02]],Plot[model,{x,xMin,xMax},PlotStyle->Red]
]]

plotFit[data,model/.fit]

Testing with another fit

sinData={{-4.8,542},{-4.7,489},{-4.5,407},{-4.4,331},{-4.2,247},{-4.0,191},{-3.9,170},{-3.7,200},{-3.5,256},{-3.4,337},{-3.2,410},{-3.1,476},{-2.9,498},{-2.7,486},{-2.6,447},{-2.4,346},{-2.3,266},{-2.1,199},{-1.9,164},{-1.8,171},{-1.6,226},{-1.5,303},{-1.3,397},{-1.1,459},{-0.97,493},{-0.81,486},{-0.65,429},{-0.48,358},{-0.32,291},{-0.16,220},{0,172},{0.16,174},{0.32,213},{0.48,280},{0.65,339},{0.81,416},{0.97,472},{1.1,477},{1.3,453},{1.5,387},{1.6,322},{1.8,248},{1.9,192},{2.1,169},{2.3,180},{2.4,228},{2.6,308},{2.7,380},{2.9,460},{3.1,506},{3.2,515},{3.4,449},{3.5,382},{3.7,280},{3.9,202},{4.0,172},{4.2,171},{4.4,209},{4.5,278},{4.7,362},{4.8,423},{5.0,465}};

The data above make up a scan from the electron edm molecular spin interferometer. The data are pairs of {applied magnetic field , population in F=1 hyperfine level}. They show the phase between two near-degenerate magnetic sublevels of the YbF molecule evolving as a small magnetic field is scanned. Fit the data to a (Sin^2) model with an offset and scale parameter in both x and y.
Hint: you'll need to provide some reasonable initial guesses to get the fit to converge.
Note: I'm sure you've got the point by now, but your solution shouldn't directly call Plot, ListPlot etc.

Solution

sinModel=a+b ((Sin[c x+d])^2);
sinFit=FindFit[sinData,sinModel,{{a,200},{b,300},{c,1.6},{d,0}},x];
plotFit[sinData,sinModel/.sinFit]
sinFit

Importing

Great. So we can fit data to models and plot the resulting fits. Next we'll go through another simple task that comes up frequently - importing and fitting data. This stuff is all rather pedestrian, but it'll come together to make something interesting in the end!

The Mathematica function for importing structured data is, perhaps unsurprisingly, named Import. We'll use it here to import some textual data and fit it.

Exercises

This command sets Mathematica's working directory to the data subdirectory of the directory that this notebook is in. [As long as you didn't move the notebook out of the unzipped archive then this should work.] [And assuming you're using Windows! If you're using another OS you might need to tweak it a little.]

SetDirectory[
DirectoryName[ToFileName["FileName"/.NotebookInformation[SelectedNotebook[]]]]<>"data\\"
];

Import, fit, plot

Import the data file "first.txt". The file has two columns of numbers delimited by whitespace. You should use the import mode "TABLE" to import this data. Write a function to fit and plot data. Use it to fit a Lorentzian (with background, amplitude, centre and width) to the data and plot the resultant fit on the same graph as the data. Return both the fit and the plot.
Hint: this might sound intimidatingly difficult, but it's nothing more that you've done already. Build the function step by step.
Hint 2: in case you don't remember, a Lorentzian is something like (1/(1+(((((x-c)/Γ)))^2))).
Hint 3: the last expression in a Module will be returned as the module's value when it's evaluated.
Hint 4: You're returning two things. Why not return them as a pair of rules {fit->someStuff,plot->thePlot}.
Hint 5: I often use strings as the first part of a rule, as it avoids clashes with previously defined symbols.

Solution

fitAndPlot[data_,model_,params_,variables_]:=Module[{f,pl},
f=FindFit[data,model,params,variables];
pl=plotFit[data,model/.f];
{"fit"->f,"plot"->pl}
]

And here it is in action:

lModel=a+b(1/(1+(((((x-c)/Γ)))^2)));
fitAndPlot[Import["first.txt","TABLE"],lModel,{a,b,{c,3},Γ},x]

Automating

Now we'll use our new functions to carry out a more complicated data analysis task. But before we get on to that, let's think a little about what we're actually doing when we write functions. The code we've written so far is pretty trivial in the grand scheme of things - we've only defined two functions - but already we can see what's happening. Effectively, by making a function that fits and plots the resulting fit, we've added that capability to Mathematica. It's almost like we're building a better Mathematica. In this next section we'll write code that automatically loads and fits hundreds of data files and looks for trends in the fit parameters. This code will use our functions defined above so we won't have to worry about how to do the fitting when we write it. We've already worried about that, solved it, and encapsulated it which will allow us to concentrate on just what's new. If we write the code for the next section carefully then we could further write code that takes hundreds of sets of hundreds of data files, and extracts trends of trends (which is starting to sound like the edm data analysis system!).

This layering up of functionality is what the arguing computer scientists I mentioned earlier call bottom-up design. It's a very powerful way to work. You extend the language layer by layer so that it becomes a more and more powerful language to solve your problem. There are a couple of key habits you need to get into to be able to program this way. First is to abstract out any emerging structure as early as you can: never be too lazy to write a function; it takes a little effort but you kind of get an exponential payback from writing a few good functions so it's almost always worth it. Second, if something feels wrong, fix it right away. For example, if you think a function should really be taking an extra parameter, or should be split into two functions then get into the habit of doing it (and modifying the rest of your code to reflect the change) before you do anything else. I find this quite a hard habit to keep, as it seems like a lot of work, but I think it's essential. You can't build a solid program on a dubious foundation - if you don't fix problems early you'll find that you'll get to a point where you're stuck and you can no longer face going back and fixing it because the problem is so deeply embedded in your code. When that happens to me I usually start writing again and resolve not to be so lazy next time!

It's quite difficult to work this way at first, especially when you find that minor problems of not knowing the language well are making things harder. But I'd urge you to keep these ideas in mind when you're writing your own code. It's worth the effort.

Right. That's enough preaching for one workbook. Back to the exercises.

Exercises

There are 376 data files in the data directory named "D1.txt" through "D376.txt". These files each contain a scan over a Lorentzian line - the data is in the same format as above. In the following exercises we'll fit all of these Lorentzians and look for a pattern in the fit parameters.

Importing and fitting a scan

First we should write a function to import and fit a scan to a Lorentzian. We've done almost all of the work already - here we'll just tie up a few ends. Your function should take the scan filename as a parameter.

Test the new function on "D75.txt". 

Hint: to keep things from getting too complicated you can assume that giving the fitter a starting parameter of 2 for the Lorentzian centre, and 0.5 for the Lorentzian width is all that is needed to make it work. You can assume that this will be true for all of the scans.

Solution

importAndFit[filename_]:=fitAndPlot[Import[filename,"TABLE"],lModel,{a,b,{c,2},{Γ,0.5}},x]
importAndFit["D75.txt"]

Generating a list of file names

Write a pure function that takes an integer, n, and returns the string "Dn.txt". Map this function over the range 1 to 376 to generate a list of filenames to analyse.
Hint: StringJoin (infix shorthand <>) and ToString might be useful here.

Solution

filenames="D"<>ToString[#]<>".txt"&/@Range[376];

Fitting the data

Now you have a function that takes a data filename and a list of filenames. Map the importing/fitting function over the list of files to generate a list of fits.
Hint: remember that a semicolon will stop Mathematica from printing out all of the results.
Note: this should take around ten to twenty seconds to run. If it takes much more or less then start to worry.

Solution

fits=importAndFit/@filenames;//AbsoluteTiming

Looking for a pattern

Now we can look for a pattern in the fit parameters. Let's assume - this example is rather contrived - that the secrets of the universe will be revealed if we plot the fitted background on the x-axis and the fitted Lorentzian amplitude on the y-axis. ListPlot these fit parameters as just described.
Hint: remember the fits are returned as lists of rules. You can use the /. operator to effect transformations with these rules.

Solution

ListPlot[{a,b}/.("fit"/.fits),AspectRatio->1]

Uncertainty estimates

In the above examples we've not considered what the uncertainties in our fitted parameters are. These uncertainties are of course essential if we are to interpret our (or others’) results correctly.

[If you don't know much about statistics then you'll probably find the following section rather hard to follow.  I'd strongly suggest spending a little time reading up on it though because it's something you need to know about. A nice, clear introduction can be found in the chapters on fitting data in Numerical Recipes. It's well worth having a decent understanding of statistics if you want to report results in your papers that are, well, you know, right. I'll try to refrain from ranting, but you might be surprised by how many papers you'll read where the authors make a hash of basic statistics. Theorists are not exempt from needing this understanding - presuming that you might want to tie your work to actual results!]

The function FindFit does not offer errors in the fit parameters, but two similar functions LinearModelFit and NonlinearModelFit (note the strange capitalisation) do. Before we investigate these functions let's quickly revise what kinds of confidence estimates we should expect to be able to get. For a linear fit (note well that linear regression applies to models that are linear in the parameters, not the variables) with errors on the data points, the errors on the fitted parameters are uniquely determined without further assumption. This is the best, and hopefully most common case. A linear fit with no errors on the points can still yield errors on the fitted parameters, but only under the assumptions that the data have similar error bars, and we've correctly chosen the form of the model underlying the data - note that this can be a very dangerous assumption. [Amazingly, this point is not understood by a number of reasonably eminent physicists I've met. Make sure that you do understand it!] For non-linear fits it gets a bit more tricky - one of the standard techniques is to find confidence ellipsoids by finding the excursions in parameter space from the optimum that change (χ^2) by a certain amount. Again this only works without further assumption if you have error bars on the data points. 

So how do Mathematica's regression functions relate to these techniques? We'll take a look at that in the following exercises and examples.

Fitting with LinearModelFit

LinearModelFit is used much like FindFit. The differences are: you give it a list of basis functions instead of a model; and it returns much more information about the fit. Let's use it to fit a straight line - I'll make up some data with slope 3, offset 5, and Gaussian y-noise with a standard deviation of 4 on the points. I'll use ErrorListPlot to plot the data with error bars.

Needs["ErrorBarPlots`"]

linData=({#1,3 #1+5+Random[NormalDistribution[0,4]]}&)/@Range[30];
ErrorListPlot[{({#1,ErrorBar[4]}&)/@linData}]

fit=LinearModelFit[linData,{1, x},x]

This FittedModel has a lot of information hidden inside it. Take a look at its help page for more details. We can extract what it calls "Standard Error"s on these parameters. Note, however, that we haven't yet let LinearModelFit know about our error bars. That means that Mathematica must be making some assumption to generate these errors. Below, we'll compare these standard errors to those that we would calculate by hand.

fit["ParameterErrors"]

Quick exercise

Make up some polynomial data with noise. Fit it with LinearModelFit.

Solution

error=20;
data=({#1,3 ((#1)^2)+5-(((#1)^3)/10)+Random[NormalDistribution[0,error]]}&)/@Range[30];
ErrorListPlot[{({#1,ErrorBar[error]}&)/@data}]
fit=LinearModelFit[data,{1,x,(x^2),(x^3)},x];

fit["BestFit"]

Linear weighted fit

To get really useful errors in the fit parameters we need to let Mathematica know what our error bars are. We do this by performing a weighted fit. The larger a point's error bar, the smaller the weight we'll give it in the fit. LinearModelFit's documentation is a bit vague as to what it will actually do, but it says that:

 " Using VarianceEstimatorFunction->(1&) and Weights->{1/((Subscript[Δy, 1])^2),1/((Subscript[Δy, 2])^2),…}, (Subscript[Δy, i]) is treated as the known uncertainty of measurement (Subscript[y, i]) and parameter standard errors are effectively computed only from the weights."
 
So let's try that. Re-running the fit with these weights (remember that I generated the data with noise on the y-axis of standard deviation 4).

fit=LinearModelFit[linData,{1,x},x,VarianceEstimatorFunction->(1&),Weights->((1/(4^2))&/@linData)];
fit["BestFit"]
fit["ParameterErrors"]

How do these estimates compare to the analytic results for the uncertainties in the fit parameters? (You can find the following formulae and their derivation in Numerical Recipes, section 15.2, available free online here).

linearFitCombs[dat_]:={s->Plus@@(1/(((#[[3]])^2)&/@dat)),sx->Plus@@((#[[1]]&/@dat)/(((#[[3]])^2)&/@dat)),sy->Plus@@((#[[2]]&/@dat)/(((#[[3]])^2)&/@dat)),sxx->Plus@@((((#[[1]])^2)&/@dat)/(((#[[3]])^2)&/@dat)),sxy->Plus@@((#[[1]] #[[2]]&/@dat)/(((#[[3]])^2)&/@dat))}

linearFit={Δ->s sxx-(sx^2),a->((sxx sy-sx sxy)/Δ),b->((s sxy -sx sy)/Δ),(Subscript[σ, a])->(Sqrt[(sxx/Δ)]),(Subscript[σ, b])->(Sqrt[(s/Δ)])};

{a,b,(Subscript[σ, a]),(Subscript[σ, b])}//.linearFit/.linearFitCombs[Join[#,{4}]&/@linData]//N

They're the same, which is reassuring! Interestingly they're similar to the values Mathematica produced without knowing about the error bars, but not the same. Once again though, I stress that Mathematica is making assumptions when deriving these errors without point error bars. We'll see how one of these assumptions could break down right now.

Different weights

What about if the data don't all have the same weight? First of all neglect the error bars, and then try fitting a line using LinearModelFit and the analytic formulae above for data that have unequally weighted points. How do the fits and errors compare?

Solution

data=({#1,-5+10 #1+Random[NormalDistribution[0,#1]]}&)/@Range[30];
ErrorListPlot[{({#1,ErrorBar[#1[[1]]]}&)/@data}]

fit=LinearModelFit[data,{1, x},x]
fit["ParameterErrors"]

fit=LinearModelFit[data,{1, x},x,VarianceEstimatorFunction->(1&),Weights->((1/((#[[1]])^2))&/@data)]
fit["ParameterErrors"]

{a,b,(Subscript[σ, a]),(Subscript[σ, b])}//.linearFit/.linearFitCombs[Join[#,{#[[1]]}]&/@data]//N

With this dataset the non-weighted results don't capture the fact that we can determine the intercept rather well Again, the weighted and "analytic" results reassuringly match up.

Playtime

Once again, here are some exercises that you can have a go at if you're feeling keen. They expand upon a few of the points we've covered above.

Do it yourself

Try doing a non-linear fit yourself by minimising (χ^2) explicitly. Solve for fit parameters that give a (χ^2)=((Subscript[χ, min])^2)+1, where ((Subscript[χ, min])^2) is the value of (χ^(2 ))at the minimum (you'll want to vary them one at a time leaving the others fixed). This should give you estimates for the errors in the fitted parameters. Try and compare these to Mathematica's parameter errors.

Solution

I'll re-use the (Sin^(2 ))interferometer data set.

chiSquare[data_,model_,parameterRules_]:=Plus@@((((model/.parameterRules)/@(#[[1]]&/@data)-(#[[2]]&/@data)))^2)

I'll numerically minimise the (χ^2) function.

sinModelFunc=a+b ((Sin[c #+d])^2)&;
sol=NMinimize[chiSquare[sinData,sinModelFunc,{}],{{a,100,200},{b,300,400},{c,1.5,1.7},{d,-1,1}}]

In the absence of errors on the data points I'll have to assume that they all have the same error and the model fits ideally. This means I should normalise my (χ^(2 ))to the number of degrees of freedom.

csNorm=((Length[sinData]-4)/(sol[[1]]));

oneSigmaCI=Solve[csNorm chiSquare[sinData,sinModelFunc,Drop[sol[[2]],{#}]]==(csNorm sol[[1]])+1,sol[[2]][[#]][[1]]]&/@{1,2}

(a/.sol[[2]])-(a/.oneSigmaCI[[1]])
(b/.sol[[2]])-(b/.oneSigmaCI[[2]])

nlFit=NonlinearModelFit[sinData,sinModel,{{a,200},{b,300},{c,1.6},{d,0}},x];
nlFit["BestFitParameters"]

nlFit["ParameterErrors"]

Bootstrapping the dataset

A very simple and general technique for estimating the uncertainty in any statistical estimator is the bootstrap (here the fitting function is the estimator - it estimates the coefficients in the fitted model). For now let's not worry too much about the details of how the bootstrap works: it's easier to do than to explain! The basic idea is to generate fake data sets from the real data set. You do this by pulling out points from the real data set at random until you have as many points in the fake data set as in the real one. You don't worry if you pull out the same point more than once. You then apply the estimator to the fake data set (in this case, fit it). Do this a few thousand times and you get a histogram of likely estimated values. Take the standard deviation of this and you've got an estimate for the error in the estimator! Easy, right?

To get you started, here's a function that generates a new data set (often called a replicate) from an original data set.

makeBootstrapReplicate[data_]:=data[[RandomInteger[{1,Length[data]},Length[data]]]]

makeBootstrapReplicate[linData]

Fitting a replicate

Write a function that generates a bootstrap replicate and fits it. The function should return the fit parameters.
Hint: because you're not interested in errors, just the fit parameters, you can use FindFit if you'd rather.

Solution

fitBootstrap[data_]:=FindFit[makeBootstrapReplicate[data],a+b x,{a,b},x]

Bootstrap estimate of parameter uncertainties

Build a table of fit results using the function you've just written. (10^4) fits should be plenty. Histogram the fitted parameter estimates. Take the standard deviation of the fitted parameter estimates to estimate the error in the fitted paramters. Compare these errors to those given above.
Hint: there's a function Histogram and a function StandardDeviation.

Solution

bsData=Table[fitBootstrap[linData],{(10^4)}];//Timing

GraphicsRow[{Histogram[a/.bsData],Histogram[b/.bsData]}]

StandardDeviation/@Transpose[({a,b}/.bsData)]

The bootstrapped results agree reasonably with both Mathematica's estimates and the analytic estimate.