`__
for more on the statistical background and the corresponding biases.)
Parameter Uncertainties and Correlations: Bootstrap and MCMC
------------------------------------------------------------
As you probably noticed, part of the output of imfit is a set of 1-sigma
parameter uncertainties for each fitted parameter in the model. These
are automatically generated when using the default (Levenberg-Marquardt)
minimizer. They’re not usually all that accurate, they assume the
uncertainties are all symmetric, and they don’t provide any information
about possible correlations or anti-correlations between different
parameter values.
If you would like a better picture of what the parameter uncertainties
and possible correlations are like, there are two options: one fast but
noisy, the other slow but detailed:
1. **Boootstrap resampling**: This involves generating a new version of
the data image by sampling from the original image with replacement
(ignoring masked pixels) and re-running the fit. Do this several
hundred (or ideally several thousand) times, and you get a
distribution of parameter values that can approximate the likelihood
(e.g., the χ2).
2. **Markov chain Monte Carlo (MCMC) analysis**: This involves computing
Markov chains consisting of sequences of sets of parameter values.
After an initial “burn-in” period, the distribution of points in
parameter space represented by a chain should converge to something
proportional to the likelihood. (The particular algorithm used by
Imfit actually runs multiple chains in parallel.)
Bootstrap Resampling Example
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
To save time, we’ll use the model *without* PSF convolution (you can of
course use PSF convolution with bootstrap resampling; it will just take
longer):
::
imfit ic3478rss_256.fits -c config_sersic_ic3478_256.dat
--mask ic3478rss_256_mask.fits --gain=4.725 --readnoise=4.3
--sky=130.14 --bootstrap 500 --save-bootstrap=bootstrap_output.dat
This will do the fit as before, print the result, and then start doing
500 rounds of bootstrap resampling and fits to the resampled data. When
it’s done (this takes about 30 seconds on a 2012 MacBook Pro with a
quad-core CPU) it will print out a summary of the best-fit parameter
values and their uncertainties; it will also save all 500 sets of
parameter values in the file bootstrap_output.dat.
This file has one column per parameter; the column names are the
parameters with numbers appended (e.g., ``X0_1``, ``n_1``) to make it
possible to distinguish different parameters when multiple versions of
the same function, or just multiple functions that have the same
parameter names, are used in the model. (I.e., all parameters for the
first function will have ``_1`` appended, all parameters from the second
will have ``_2`` appended, etc.)
In the ``python/`` subdirectory of the main Imfit package there are a
couple of Python modules: imfit_funcs.py and imfit.py. The latter has a
simple function to read in the bootstrap-resampling output file
(``imfit.GetBootstrapOutput``), which will return a list of parameter
names and a 2D Numpy array with the full set of parameter values.
There are many possible ways of analyzing the bootstrap-resampling
output. One thing you can do, if the model is not *too* complicated, is
make a scatterplot matrix (a.k.a. corner plot) of the parameters. The
Python package `corner.py `__
can be used for this; here’s a quick-and-dirty example that also uses
the ``imfit.GetBootstrapOutput`` function:
::
>>> import imfit, corner
>>> columnNames, bootstrapResults =
imfit.GetBootstrapOutput("bootstrap_output.dat")
>>> corner.corner(bootstrapResults, labels=columnNames)
The result is shown below.
.. raw:: html
.. image:: ./fig4_for_tutorial.png
**Figure 4:** Scatterplot matrix of parameter values from 500 rounds of
bootstrap resampling fits to the IC 3478 *r*-band image (Sérsic model,
no PSF convolution). Note the clear correlations between the Sérsic
model parameters (n, r_e, I_e).
MCMC Example
~~~~~~~~~~~~
MCMC analysis uses a separate program called ``imfit-mcmc``. You can run
it with the following command (note that it’s identical to the regular
``imfit`` command, except for the option that specifies the root name
for output files):
::
imfit-mcmc ic3478rss_256.fits -c config_sersic_ic3478_256.dat
--mask ic3478rss_256_mask.fits --gain=4.725 --readnoise=4.3
--sky=130.14 --output=mcmc_ic3478r
**Warning:** this will take several minutes! (On my 2012 MacBook Pro
with a quad-core Intel i7 CPU, it takes about eight or ten minutes.)
Various updates will be printed as the program runs. Once a trial
“burn-in” phase is over, ``imfit-mcmc`` will test for possible
convergence of the chains every 5,000 generations by looking at the last
half of each chain. If convergence is detected, the program will quit;
otherwise, it will quit when it reaches 100,000 generations. (These
values can be changed with command-line options.)
When it’s done, you will have *seven* output text files, named
mcmc_ic3478r.1.txt, mcmc_ic3478r.2.txt, etc., one for each of the
individual chains. (By default, the total number of chains is equal to
the number of free parameters in the model.) Each is similar to the
bootstrap-resampling output file in format, with one column for each
parameter in the model (plus some extra bookkeeping columns that you can
ignore unless you’re interested in details of the MCMC process), and one
row for each generation in the chain; each chain will have several tens
of thousands of generations.
The ideal thing to do is probably to take the last half of each chain
and combine them all into one gigantic set of parameter values. There’s
a Python function for that in python/imfit.py, which returns the same
kinds of output as imfit.GetBootstrap (i.e., a list of parameter names
and a 2D Numpy array). Here’s an example of using that, and then making
a scatterplot matrix with the corner.py module, just as we did for the
bootstrap output:
::
>>> import imfit, corner
>>> columnNames, allchains = imfit.MergeChains("mcmc_ic3478r",
secondHalf=True)
>>> corner.corner(allchains, labels=columnNames)
The result is shown below.
.. raw:: html
.. raw:: html
.. image:: ./fig5_for_tutorial.png
**Figure 5:** Scatterplot matrix of parameter values from Markov chain
Monte Carlo analysis of the IC 3478 *r*-band image (Sérsic model, no PSF
convolution). Note the strong correlations between the Sérsic model
parameters (n, r_e, I_e), and the weaker correlation between r_e and
ellipticity and between X0 and Y0. Since this plot is based on about
300,000 samples, it is considerably less noisy than the version based on
500 rounds of bootstrap resampling in `Figure 4 <#fig4>`__.
.. raw:: html