PHYCS 6730 Lab Exercise: Nonlinear Optimization: Chi Square

These exercises experiment with the problems of chi square fits to nonlinear data. The difficulties here are common to general nonlinear optimization problems. The files you need are in ~p6730/exercises/nonlinear_chisq.

Exercise 1 The spectral data

Using your favorite plotting utility, plot the files peak1 and peak3. The data are simulations of Gaussian spectral lines with random fluctuations added to simulate noise. Notice that there are three numbers on each line. They are x, y, and the standard deviation in y. So it helps to plot the points with error bars. These are artificially generated data intended to represent a Gaussian peak and two superimposed Gaussian peaks, respectively, each with artificially added noise. It is difficult to tell the two peaks in peak3.

In these exercises we will be exploring the chi square surface as a function of fit parameters, where the fitting function is (Maple notation)

   yexp(w) = a1 + a2*exp(-1/2*(w - a3)^2/a4^2) + a5*exp(-1/2*(w - a6)^2/a7^2)

That is, the fitting function has a baseline constant a1 and a sum of Gaussians. Each of the Gaussian peaks is specified by three parameters: the amplitudes (a2 and a5) the peak positions (a3 and a6) and the widths (a4 and a7). If we have only a single peak, we drop the second term and keep only four parameters.

Nothing to hand in for this exercise.

Exercise 2. Single peak. Maple contour plot

In the file mpeak1 I have preconstructed a Maple session for your viewing pleasure. Your task is to run it, play with it, and understand the results.

Start xmaple and load the file using

   read mpeak1;

(Don't forget the semicolon!)

The best fit values for the parameters are defined. They are close to the values used to create the data: (0.2, 1, 140, 8) before we started adding random noise to the values.

Find the formula for chi square in the Maple commands you loaded. What parameters does it depend on?
Evaluate chi square at the best fit values. You should get 36.5 for this case. For this purpose, notice that chisq is defined to be a Maple function (procedure) that takes four variables.
Evaluate chi square at the values used to create the data. (See above.) You should get a rather higher value. That's life! Does the higher chi square value exceed the number of points?
Since we are dealing with a four-dimensional space, we have to view the chisquare surface in slices. So let's draw the chisquare contours in the two linear parameters a1 and a2 with the last two parameters fixed at their best fit values. Use the Maple command
```
with(plots);
contourplot(chisq(a1,a2,a3best,a4best),a1=0.20..0.24,a2=1.1..1.4,contours=[37.,38.,39.,40.,41.]);
```
(If you have trouble copying and pasting from your browser to Maple, see the file commands1. It contains a text version of the command that is easier to copy.)
See if you can locate the minimum and compare with the best fit values. Notice how this command specifies a plotting range as well as contour levels. You are encouraged to play with these settings here and throughout this lab exercise.
Now let's see how things change when the dependence on the parameter is nonlinear. We fix a1 and a3 at the best fit values and vary the other two parameters.
```
contourplot(chisq(a1best,a2,a3best,a4),a2=0..3,a4=1..20,contours=[40.,100.,150.,200.,400.]);
```
Here we have specified a coarser set of contours. What qualitative difference do you see in the shapes of the contours?
Next let's look at the a3, a4 plane with the first two parameters fixed at their best fit values.
```
contourplot(chisq(a1best,a2best,a3,a4),a3=100..150,a4=1..20,contours=[40.,100.,150.,200.,400.]);
```
Then let's zoom in on the minimum.
```
contourplot(chisq(a1best,a2best,a3,a4),a3=137..141,a4=5.5..7.0,contours=[37.,38.,39.,40.,41.]);
```
Why should the contours look more elliptical when we are close to the minimum?

Exercise 3. Double peak - double jeopardy

Here we play with a seven parameter space. Reset your Maple session and load the mpeak3 set:

reload;
read mpeak3;

These data were generated with the parameter values (0.2,1,140,8,0.6,147,6). Unfortunately, the noise is enough to obscure the peaks and the best values turned out to be (0.194,0.099,125,5.6,1.6,143,7.5). So the best fit really went off base here.

Note that the values of a1gen, etc are preset to their intended values, instead of the rather strange best fit values.

Evaluate chi square at the intended and best fit points.
As before, let's examine the dependence on linear parameters a2 and a5 in the vicinity of the intended parameter values.
```
contourplot(chisq(a1gen,a2,a3gen,a4gen,a5,a6gen,a7gen),a2=0..2,a5=0..2,contours=[35,38,41,44]);
```
Is there a well-defined minimum close to the intended values of a2 and a5?
Now let's see what happens to the contours if we force a3 = a6 and a4 = a7.
```
contourplot(chisq(a1gen,a2,142,8,a5,142,8),a2=-1..2,a5=0..2,contours=[35,38,41,44]);
```
How would you determine the minimum in this case? Can you explain why the contours look like straight lines? Consider how the fitting function depends on the two parameters in this case.
Now let's see how the minimization of chi square got lost. We try forcing the location of the first peak to be a3 = 124, close to the best fit value for this parameter. See how the contours in the a2, a5 plane change in this case:
```
contourplot(chisq(a1gen,a2,124,8,a5,142,8),a2=-1..2,a5=0..2,contours=[35,38,41,44]);
```
Where is the minimum now?
So let's set parameter a5 to the new best value and explore the a2, a3 plane. In this plane we are allowing the first peak position and amplitude to vary, while keeping the other parameters fixed. Notice that the parameters a5, a6, a7 are carrying the entire peak in the data and we are now asking whether there is any place for a second peak with a width 8.
```
contourplot(chisq(a1gen,a2,a3,8,1.4,142,8),a2=-1..2,a3=100..300,contours=[35,38,41,44]);
```
Where are the minima now? Can you explain how the minimization process got lost? What are the characteristics of the data that led to this muddle?