Data

1   # number of galaxies
# x y err
0.000000e+00 0.000000e+00 3.000000e-03
# the next three lines are vestigial
0.0 10000.0
0.0 10000.0
0.0 10000.0

Next come lines to specify the image properties. The code can handle different sets of images corresponding to different sources. You first specify the number of distinct sources. Then for each source you specify the number of images followed by the properties of the images.

1   # number of sources

4   # number of images of this source
# x y flux err(x) err(flux) tdel err(tdel)
2.175513e-01 -6.543659e-01 -1.139633e+00 3.000000e-03 5.698163e-02 6.481701e-01 0.000000e+00
2.669991e-01 1.063528e+00 -8.091063e+00 3.000000e-03 4.045531e-01 1.049748e+00 0.000000e+00
6.794387e-01 9.189104e-01 8.494827e+00 3.000000e-03 4.247413e-01 1.185709e+00 0.000000e+00
-1.260121e+00 2.319804e-01 3.003357e+00 3.000000e-03 1.501678e-01 1.408355e+00 0.000000e+00

The data file is loaded with the data command.

Some comments:

• Fluxes: the units are arbitrary; the code will scale the source flux appropriately to find the best fit.
• Parities: if you know (or can guess) the parities, you can give the flux an appropriate sign.
• Time delays: any time delay whose errorbar is 0 is ignored.
• Error ellipses: the code can handle position uncertainties represented by error ellipses; just replace the position errorbar by three numbers giving the semi-major axis, semi-minor axis, and position angle of the error ellipse, and then load the data file with the data2 command.
• Central images: if central images have not been observed, you can omit them here and tell the code to ignore predicted central images using omitcore (see below).

Goodness of fit

The code tracks and reports the separate contributions from the image positions, image fluxes, galaxy position, and time delays (in addition to the total).

The code has two options for evaluating the chisq contribution from the image positions.

• image plane chisq: The code uses the tiling algorithm to solve the lens equation and find the model images. It then compares the observed and predicted positions in the image plane.
• source plane chisq: The code maps each observed image back to the source plane, and computes the scatter among the source positions. It includes a factor of the magnification tensor so the source plane chisq approximates the image plane chisq.

However, there is one circumstance in which the image plane chisq is required: when the number of images is an issue. The image plane chisq checks that the number of model images equals the number of observed images; but the source plane chisq does not. When modeling standard double or quad lenses, this is rarely a problem, because if a model fails to predict the right number of images it will probably be bad for other reasons. When modeling unusual lenses, by contrast, it may be important to ensure that the model predicts the correct number of images. (PMN J0134 provides a nice example.)

Minimizing the chisq

The code offers two routines from Numerical Recipes for finding minima of a nonlinear function in a multi-dimensional parameter space:

• amoeba: downhill simplex method
• powell

Both amoeba and powell involve keeping track of N+1 points in the N-dimensional parameter space. You can specify multiple starting points if you like (in the startup or setlens command). If you specify fewer than N+1 models, the code will fill in the rest automatically.

Non-linear least squares fitting always faces the danger of getting stuck in local minima, or hung up in long, narrow valleys. My strategy for handling this problem is to wait until the optimization routine has stopped, and then restart it. The code keeps the best model so far, but generates the other N starting points from scratch. You can tell the code to do this using the restart variable. In my experience, restarting 3-5 times is sufficient for most problems.

I should say that the code actually does two nested optimizations:

• 1. With the lens model fixed, the code optimizes the source position and flux.
• 2. Then optimize the lens model (doing step #1 for each point in the model parameter space).

Exploring parameter space (I): Random search

When you first start to model a new lens, you probably have no idea what parameter values are reasonable. You need to make a first survey of parameter space to figure out where the reasonable models are. I have tried various strategies over the years, and my favorite approach now is to pick random points in parameter space as starting points for the optimization. Each realization will lead to some local minimum; if there is a single, global minimum, all the realizations will lead there. Working with the source plane chisq, this process is very fast. The following input file gives an example.

# ignore images within 0.01 of a galaxy
set omitcore = 0.01
# read data
data mock3-q1.lens
# use source plane chisq, so we do not need grid
set chimode = 0
set gridflag = 0
set checkparity = 0
# to handle local minima and narrow valleys
set restart = 3
# specify the model class
setlens 1 1
  alpha 1 0 0 0.3 10 0.1 10 0 0 1
  1 1 1 1 1 1 1 0 0 0
# do 10 optimizations starting from random points
randomize 10 sieg-ran
  0.5 1.5
  -0.01 0.01
  -0.01 0.01
  0.0 0.7
  -90.0 90.0
  0.0 0.3
  -90.0 90.0
# take the best model found so far and reoptimize
setlens sieg-ran.start
optimize sieg
#
quit

Here is one way to look at the results. The open squares show the random starting points in the plane of ellipticity and its position angle (left), and shear and its PA (right). The red and green colors show two of the 10 independent realizations. We see that both realizations converged to the same model, shown by the blue point.

Here is the final model:

LENS PARMS:
  alpha 1.000019e+00 -1.053656e-05 -1.575588e-06 2.998274e-01 -7.860810e-03 1.000592e-01 2.499274e+01 0.0 0.0 1.000000e+00
SOURCE PARMS:
  ptsrc 9.998177e-01 -1.288372e-01 1.867252e-01 0.0 0.0 0.0 0.0 0.0
CHISQ: 7.500291e-05 5.795483e-05 4.425894e-06 0.000000e+00 1.262219e-05 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 # tot pos flux tdel gal plim crv ring weak

Source #1
  posn: -1.288372e-01 1.867252e-01
  flux: 9.998177e-01
  h: 7.000000e-01
  tscale: 1.000000e+00
images:
0.2176 -0.6544 (3.0e-03) -1.1396 ( 5.70e-02) 0.65 (0.0e+00) -> 0.2175 -0.6544 -1.1396 0.7613
0.2670 1.0635 (3.0e-03) -8.0911 ( 4.05e-01) 1.05 (0.0e+00) -> 0.2670 1.0635 -8.0917 0.2577
0.6794 0.9189 (3.0e-03) 8.4948 ( 4.25e-01) 1.19 (0.0e+00) -> 0.6794 0.9189 8.4949 0.2520
-1.2601 0.2320 (3.0e-03) 3.0034 ( 1.50e-01) 1.41 (0.0e+00) -> -1.2601 0.2320 3.0032 0.0000

In this case I modeled a mock lens created from an SIE galaxy with ellipticity and shear, so the SIE+shear model fits the data and recovers the input parameters perfectly.

Exploring parameter space (II): Controlled searches

Always explore parameter space! You MUST understand the uncertainties in your model results.

Lensmodel provides a variety of ways to take slices through parameter space. First, you can plot chisq versus a single parameter, allowing some or all of the other parameters to vary:

varyone igal ip lo hi #steps outbase

You specify the parameter to explore by giving the galaxy index igal and the parameter number ip. Then you specify the range of parameter values and number of steps to take. Here are examples:

varyone 1 1 0.98 1.02 21 vary-b
varyone 1 4 0.20 0.40 21 vary-e

Here are the results plots of chisq versus the lensing strength (left) and ellipticity (right).

You can also take 2-d slices as follows:

varytwo   igal ip lo hi #steps   jgal jp lo hi #steps   outbase

You specify the indices and ranges of the two parameters, Here is an example with ellipticity and shear.

varytwo   1 4 0.20 0.40 31   1 6 0.05 0.15 31   vary-eg

And here is a plot of the 1, 2, and 3-sigma contours.

Exploring parameter space (III): Markov chain monte carlos

MCMC is a powerful and increasingly popular method for characterizing the posterior probability distribution for the model. The basic idea is to find the best-fit model, and then take random steps around that model. If you do the random search just right, you wind up sampling from the posterior probability distribution, no matter how complicated it is, and no matter whether you can actually characterize it in any simple way. Now is not the time to go into theoretical details; let me just illustrate how it can be done in lensmodel:

stepsize
  0.003 0.003 0.003 0.01 0.1 0.01 0.1 0 0 0
mcmc sieg-mcmc 10

The stepsize line specifies an initial guess for the range of each parameter. In the mcmc line, the second argument is the number of chains to run.

Here are sample results for the lens galaxy position. The jagged lines show five of the random chains, while the points represent the final sample of points drawn from the posterior distribution.

Here for the ellipticity and its position angle (left), and the shear and its position angle (right).

Or if you prefer to see the ellipticity and shear plotted against each other (left), and the two position angles (right).

In case you are curious how the MCMC results compare with the chisq contours from varytwo, here are the results in the plane of ellipticity and shear.