Fitting Models to Data

As we have seen, we can use HoloPy to perform Scattering Calculations from many types of objects. Here, the goal is to compare these calculated holograms to a recorded experimental hologram, and adjust the parameters of the simulated scatterer to get a good fit for the real hologram.

A Simple Least Squares Fit

We start by loading and processing data using many of the functions outlined in the tutorial on Loading Data.

import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path, load_average
from holopy.core.process import bg_correct, subimage, normalize
from holopy.scattering import Sphere, Spheres, calc_holo
from holopy.inference import prior, ExactModel, CmaStrategy, EmceeStrategy

# load an image
imagepath = get_example_data_path('image01.jpg')
raw_holo = hp.load_image(imagepath, spacing = 0.0851, medium_index = 1.33,
                         illum_wavelen = 0.66, illum_polarization = (1,0))
bgpath = get_example_data_path(['bg01.jpg','bg02.jpg','bg03.jpg'])
bg = load_average(bgpath, refimg = raw_holo)
data_holo = bg_correct(raw_holo, bg)

# process the image
data_holo = subimage(data_holo, [250,250], 200)
data_holo = normalize(data_holo)

Next we define a scatterer that we wish to model as our initial guess. We can calculate the hologram that it would produce if it were placed in our experimental setup, as in the previous tutorial on Scattering Calculations. Fitting works best if your initial guess is close to the correct result. You can find guesses for x and y coordinates with center_find(), and a guess for z with propagate().

guess_sphere = Sphere(n=1.58, r=0.5, center=[24,22,15])
initial_guess = calc_holo(data_holo, guess_sphere)
hp.show(data_holo)
hp.show(initial_guess)

Finally, we can adjust the parameters of the sphere in order to get a good fit to the data. Here we adjust the center coordinates (x, y, z) of the sphere and its radius, but hold its refractive index fixed. By default fit() will adjust all parameters, so we can omit the argument if that is what we want.

parameters_to_fit = ['x', 'y', 'z', 'r']
fit_result = hp.fit(data_holo, guess_sphere, parameters=parameters_to_fit)

The fit() function automatically runs calc_holo() on many different sets of parameter values to find the combination that gives the best match to the experimental data_holo. We get back a FitResult object that knows how to summarize the results of the fitting calculation in various ways, and can be saved to a file with hp.save`() :

best_fit_values = fit_result.parameters
initial_guess_values = fit_result.guess_parameters
best_fit_sphere = fit_result.scatterer
best_fit_hologram = fit_result.hologram
best_fit_lnprob = fit_result.max_lnprob
hp.save('results_file.h5', fit_result)

If we look at best_fit_values or best_fit_sphere, we see that our initial guess of the sphere’s position of (24, 22, 15) was corrected to (24.16, 21.84, 16.35). Note that we have achieved sub-pixel position resolution!

Customizing the model

Sometimes you might want a bit more control over how the parameters are varied. You can customize the parameters with a Model object that describes parameters as Prior objects instead of simply passing in your best guess scatterer and the names of the parameters you wish to vary. For example, we can set bounds on the coordinate parameters and use a Gaussian prior for the radius - here, with a mean of 0.5 and standard deviation of 0.05 micrometers.

x = prior.Uniform(lower_bound=15, upper_bound=30, guess=24)
y = prior.Uniform(15, 30, 22)
z = prior.Uniform(10, 20)
par_sphere = Sphere(n=1.58, r=prior.Gaussian(0.5, 0.05), center=[x, y, z])
model = ExactModel(scatterer=par_sphere, calc_func=calc_holo)
fit_result = hp.fit(data_holo, model)

Here we have used an ExactModel which takes a function calc_func to apply on the Scatterer (we have used calc_holo() here). The ExactModel isn’t actually the default when we call fit() directly. Instead, HoloPy uses an AlphaModel, which includes an additional fitting parameter to control the hologram contrast intensity - the same as calling calc_holo() with a scaling argument. You can fit for the extra parameters in these models by defining them as Prior objects. Likewise, if the scattering theory you are using requires fittable parameters (such as the lens_angle for the MieLens theory or the spherical_aberration for the AberratedMieLens theory), you can fit for these by defining them as Prior objects as well.

The model in our example has read in some metadata from data_holo (illumination wavelength & polarization, medium refractive index, and image noise level). If we want to override those values, or if we loaded an image without specifying metadata, we can pass them directly into the Model object by using keywords when defining it.

Advanced Parameter Specification

You can use the Model framework to more finely control parameters, such as specifying a complex refractive index :

n = prior.ComplexPrior(real=prior.Gaussian(1.58, 0.02), imag=1e-4)

When this refractive index is used to define a Sphere, fit() will fit to the real part of index of refraction while holding the imaginary part fixed. You could fit it as well by specifying a Prior for imag.

You may desire to fit holograms with tied parameters, in which several physical quantities that could be varied independently are constrained to have the same (but non-constant) value. A common example involves fitting a model to a multi-particle hologram in which all of the particles are constrained to have the same refractive index, but the index is determined by the fitter. This may be done by defining a parameter and using it in multiple places. Other tools for handling tied parameters are described in the user guide on The HoloPy Scatterer.

n1 = prior.Gaussian(1.58, 0.02)
sphere_cluster = Spheres([
Sphere(n = n1, r = 0.5, center = [10., 10., 20.]),
Sphere(n = n1, r = 0.5, center = [9., 11., 21.])])

Transforming Priors

Sometimes you might want to apply mathematical operations to transform one prior into another, for example in the following use cases:

  • You want two parameters to vary together with values that are related but unequal, such as the length and radius of a cylinder with known aspect ratio, or the z-coordinates of two vertically stacked particles.
  • You want a parameter that is distrbuted according to some other distribution, such as a cylinder axis evenly distributed in spherical coordinates or a sphere with uniformly distributed volume (not radius).
  • You want to reparamaterize a problem to reduce covariances between fitting parameters, such as finding the separation distance between two closely interacting particles or to find the positions of particles confined to a plane with a slight tilt as compared to the x-y plane.

If you have an explicit transformation function you can use it to define a TransformedPrior object, for example to define a polar angle with the appropriate distribution we might do:

def uniform2polar(u):
    # we want polar angle \theta distributed according to sin(\theta)
    # Use inverse transform sampling, correct for \theta in [0, pi]
    symmetric_polar = np.arcsin(u)
    return symmetric_polar + np.pi/2
director = prior.Uniform(-1, 1, name='director')
polar_angle = prior.TransformedPrior(uniform2polar, director, name='polar')

You can use TransformedPrior objects when defining a Scatterer just like regular priors. They share some attributes with basic priors as well, such as the TransformedPrior.sample() method, but you cannot directly calculate probabilities or log-probabilities of a TransformedPrior taking on a particular value.

Besides explicitly defining them, you can also create TransformedPrior objects by using numpy ufuncs and built-in operators on priors:

sphere_area = prior.Uniform(1, 2, name='base_area')
diameter = np.sqrt(sphere_area / np.pi)
radius = diameter / 2
shell = radius + 0.1

All of our derived objects are TransformedPrior objects, even though we didn’t explicitly define them that way. If we use more than one of sphere_area, diameter, radius, or shell in a fitting or inference calculation, they will all be derived from a single parameter (sphere_area in this case) even though they take on different values. Note that we could have expressed shell in one line if we didn’t care about the intermediate values:

shell = np.sqrt(prior.Uniform(1, 2, name='base_area')) / 2 + 0.1
shell.name = 'shell'

It’s always a good idea to assign your priors names when working with TransformedPrior objects to keep track of their relationships. Parameter names will be generated if none are provided but they might not be very informative. To help with naming, you can even assign names to TransformedPrior objects when using numpy ufuncs!

diameter = np.sqrt(sphere_area / np.pi, name='diameter')

Bayesian Parameter Estimation

Often, we aren’t just interested in the best-fit (MAP) parameter values, but in the full range of parameter values that provide a reasonable fit to an observed hologram. This is best expressed as a Bayesian posterior distribution, which we can sample with a Markov Chain Monte Carlo (MCMC) algorithm. The approach and formalism used by HoloPy are described in more detail in [Dimiduk2016]. For more information on Bayesian inference in general, see [Gregory2010].

A sampling calculation uses the same model and data as the fitting calculation in the preceding section, but we replace the function fit() with sample() instead. Note that this calculation without further modifications might take an unreasonably long time! There are some tips on how to speed up the calculation further down on this page.

The sample() calculation returns a SamplingResult object, which is similar to the FitResult returned by fit(), but with some additional features. We can access the sampled parameter values and calculated log-probabilities with SamplingResult.samples and SamplingResult.lnprobs, respectively. Usually, the MCMC samples will take some steps to converge or “burn-in” to a stationary distribution from your initial guess. This is most easily seen in the values of SamplingResult.lnprobs, which will rise at first and then fluctuate around a stationary value after having burned in. You can remove the early samples with the built-in method SamplingResult.burn_in(), which returns a new SamplingResult with only the burned-in samples.

Customizing the algorithm

The fit() and sample() functions follow algorithms that determine which sets of parameter values to simulate and compare to the experimental data. You can specify a different algorithm by passing a strategy keyword into either function. Options for fit() currently include the default Levenberg-Marquardt (strategy="nmpfit"), as well as cma-es (strategy="cma") and scipy least squares (strategy="scipy lsq"). Options for sample() include the default without tempering (strategy="emcee"), tempering by changing the number of pixels evaluated (strategy="subset tempering"), or parallel tempered MCMC (strategy="parallel tempering") [not currently implemented]. You can see the available strategies in your version of HoloPy by calling hp.inference.available_fit_strategies or hp.inference.available_sampling_strategies.

Each of these algorithms runs with a set of default values, but these may need to be adjusted for your particular situation. For example, you may want to set a random seed, control parallel computations, customize an initial guess, or specify hyperparameters of the algorithm. To use non-default settings, you must define a Strategy object for the algorithm you would like to use. You can save the strategy to a file for use in future calculations or modify it during an interactive session.

cma_fit_strategy = CmaStrategy(popsize=15, parallel=None)
cma_fit_strategy.seed = 1234
hp.save('cma_strategy_file.h5', cma_fit_strategy)
strategy_result = fit(data_holo, model, cma_fit_strategy)

In the example above, we have adjusted the popsize hyperparameter of the cma-es algorithm, prevented the calculation from running as a parallel computation, and set a random seed for reproducibility. The calculation returns a FitResult object, just like a direct call to fit().

Similarly, we can customize a MCMC computation to sample a posterior by calling sample() with a EmceeStrategy object. Here we perform a MCMC calculation that uses only 500 pixels from the image and runs 50 walkers each for 2000 samples. We set the initial walker distribution to be one tenth of the prior width. In general, the burn-in time for a MCMC calculation will be reduced if you provide an initial guess position and width that is as close as possible to the eventual posterior distribution. You can use Model.generate_guess() to generate an initial sampling to pass in as an initial guess to your EmceeStrategy object.

nwalkers = 50
initial_guess = model.generate_guess(nwalkers, scaling=0.1)
emcee_strategy = EmceeStrategy(npixels=500, nwalkers=nwalkers,
    nsamples=2000, walker_initial_pos=initial_guess)
hp.save('emcee_strategy_file.h5', emcee_strategy)
emcee_result = hp.sample(data_holo, model, emcee_strategy)

Random Subset Fitting

In the most recent example, we evaluated the holograms at the locations of only 500 pixels in the experimental image. This is because a hologram usually contains far more information than is needed to estimate your parameters of interest. You can often get a significantly faster fit with little or no loss in accuracy by fitting to only a random fraction of the pixels in a hologram.

You will want to do some testing to make sure that you still get acceptable answers with your data, but our investigations have shown that you can frequently use random fractions of 0.1 or 0.01 with little effect on your results and gain a speedup of 10x or greater.