Overview of fitter
My fitter is divided into two primary parts and several auxiliary programs.
Each will be described here.
Marquardt-Levenberg Fitter with Bayes Priors
This program is based on the Levenberg-Marquardt method [1].
Briefly, this is a method to vary smoothly between the steepest descent method (used far from the minimum) and the inverse-Hessian method (used near the minimum).
One downside of the method is that it requires the derivative of the fit function which, as in our case, may be difficult to compute.
A rough outline of the program is as follows:
Driver script
The fitting program is extremely versatile, and therefore requires a large number of input values.
Typing these in by hand (or even copy-pasting them in) quickly becomes tedious.
It's also often useful to try varying a specific parameter over a range of values.
To ease the burden, I've written a script that calls the main fitter.
Here's a brief list of its capabilities:
- Ability to loop over various quantities including:
- momenta (000, 100, 110, etc)
- tmin, tmax
- bayes (y,n)
- data SVD cut
- prior widths
- datasets
- from minimum correlator to a range of maximum correlators (typically minimum is 0 and maximum range is 0,2 so it fits to local only, then local+1S, then local+2S)
- number of energy levels
- propagators (pi, ro, b1, a0, a1, hlpi, hlro)
- A good memory:
- Starting guesses for priors are read from a database file
- Intermediate results are fed back to the fitter as starting guesses
- Ability to summarize results:
- If (at least) the pi, ro, and b1 propagators are included in the fit, it will generate a charmonium spectrum mass splitting plot (generated from M1 values).
-
Detailed Descriptions
Here is a more complete description of each stage of the process:
- Read user inputs
- Which model function? Choices are:
- MultiState
- MultiState_Diag
- MultiState_ExpCosh
- MultiState_ExpCosh_Diag
- MultiState_ExpCosh_Diag_Log
- MultiOsc_ExpCosh_Diag_Log
- MultiState_ExpCosh_ELog
- Input correlators. For each correlator, the following information is requested:
- filename
- fraction of configurations to use in analysis
- fit range (tmin and tmax). If tmax <= half the timeslices then the data is folded.
- regular statistics or jackknife statistics
- source/sink values (ie, "0 0" for local-local, "1 0" for 1S-local)
- Value for rho (off-diagonal damping parameter) for the data-correlation matrix
- SVD cutoff for the data correlation matrix inversion
- Are we doing Bayes statistics? (y/n)
- How many energy levels to fit to (osc/non-osc are considered independent)?
- For each energy level, guess for energy (and prior, if we're doing Bayes)?
- How many correlators are we fitting to?
- We can automagically guess the coefficients, or input them manually. (The guesses are based on the propagator value at timeslice 1.) If we input manually, we ask for the guess (and prior, if doing Bayes)
- SVD cut for alpha matrix (Marquardt method) inversion
- Number of bootstraps to perform (0 for none)
- Method for varying energy/coefficient priors when bootstrapping (Full variation, 1-sigma maximum variation, or No variation)
- Filename to store bootstrap data
- Perform the fit
At this point, all the data is stored in a variable called alldata. The model function gets stored in model, and the function we want to minimize is stored in chi2func. The minimizer is always set up as minimizer = Levenberg_Marquardt(model, chi2func, svd_cut_ai). We perform the minimization by executing params, alpha_inv, chi2 = minimizer.full(), which performs the following steps:
- Initializes variables, eg. starting chi2, lamda (intentional misspelling due to lambda being a reserved word in python), alpha
- Attempts a minimization step (maximum of 200 steps)
- If step is good, check to see if chi2 was reduced by 1e-9 (absolute) or 1e-6 (relative). If so, then break (unless this is our first step). Otherwise, decrease lamda by factor of 2 and repeat.
- If step is bad, increase lamda by factor of 2 and try again.
- After minimizer converges or fails (200 steps without converging), restart it 10 times (to make sure we've found the minimum). These restarts shouldn't cost much since we should already be close to the minimum.
- Plot the resulting fit (currently broken)
If this were working, it would display a plot showing the fit curve going through the data. I'll get this working with the leading exponential subtraced off when I have time.
- Print the results
Results are printed in the form:
chi^2/dof = 180.894 / 32 = 5.653 (2.397 from Bayes priors)
parameter 0 is 1.53867e+00 +- 3.95114e-04 (0%) prior was 1.54000e+00 +- 5.00000e-02 (chi2+=0.001)
parameter 1 is -2.41115e-01 +- 7.03653e-03 (2%) prior was -6.73345e-01 +- 3.82743e-01 (chi2+=1.275)
parameter 2 is 3.25720e-01 +- 4.64019e-03 (1%) prior was 0.00000e+00 +- 1.00000e+00 (chi2+=0.106)
parameter 3 is 1.00731e+00 +- 2.46785e-03 (0%) prior was 0.00000e+00 +- 1.00000e+00 (chi2+=1.015)
- If bootstraps have been requested:
- Else:
References
1 Press et. al., Numerical Recipes, Section 15.5 "Nonlinear Models"