Bayesian inference of Parameter Values for Multiple Spheres

[rbnvrw]

First of all, thank you for making this code publicly accessible and thank you for all the work that has gone into this.

I was following the tutorial on "Bayesian inference of Parameter Values" and I am trying to adapt it to multiple spheres. MWE:

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, calc_holo
from holopy.inference import prior, AlphaModel, tempered_sample
from holopy.scattering import Spheres

s1 = Sphere(center=(5, 3, 5), n = 1.59, r = 1)
s2 = Sphere(center=(5, 5, 5), n = 1.59, r = 1)
s3 = Sphere(center=(5, 7, 5), n = 1.59, r = 1)
collection = Spheres([s1, s2, s3])

medium_index = 1.33
illum_wavelen = 0.405
illum_polarization = (1,0)
detector = hp.detector_grid(shape = 200, spacing = .05)
holo_405 = calc_holo(detector, collection, medium_index, illum_wavelen, illum_polarization)

# Set up the prior
pg5 = prior.Gaussian(4.97, .08)
pn = prior.Gaussian(1.55, .1)
pr = prior.BoundedGaussian(0.98, .05, 0, 5)
s1 = Sphere(center=(pg5, prior.Gaussian(2.97, .08), pg5), n = pn, r = pr)
s2 = Sphere(center=(pg5, pg5, pg5), n = pn, r = pr)
s3 = Sphere(center=(pg5, prior.Gaussian(6.97, .08), pg5), n = pn, r = pr)
priors = Spheres([s1, s2, s3])

noise_sd = holo_405.std()
model = AlphaModel(priors, noise_sd=noise_sd, alpha=1)
result = tempered_sample(model, holo_405, nwalkers=10, samples=100, max_pixels=200*200*0.0005)

result.values()

If I run this code, I get a lot of OverlapWarnings and the following output:

{'0:Sphere.center[1]': UncertainValue(value=2.991103114359165, plus=0.27333182893886043, minus=0.14208225419107867, n_sigma=1),
 '2:Sphere.center[1]': UncertainValue(value=6.8950442696398495, plus=0.1741085966110747, minus=-0.007508585247170352, n_sigma=1),
 'Sphere.n': UncertainValue(value=1.6713015019771655, plus=0.13013135982273671, minus=0.1576104705835486, n_sigma=1),
 'Sphere.r': UncertainValue(value=0.9487928083946258, plus=0.341028184657435, minus=0.0608332119720173, n_sigma=1),
 ']': UncertainValue(value=5.031485549049074, plus=0.10400240255691262, minus=0.09450660658172616, n_sigma=1)}

So, it seems the code is able to correctly determine the y coordinates, but for the other parameters I see no values. I could not find anything in the documentation about multiple spheres, however in your paper ("A Bayesian approach to analyzing holograms of colloidal particles", Dimiduk & Manoharan, 2016) you write in the conclusions "In the future, we aim to apply this technique to more complex scattering models such as the superposition solution for multiple spheres". Is this supported yet?

[tdimiduk]

You have accidentally wandered into some behavior that is not as well documented as it should be.

If you use the same prior object (pg5 in your case) in multiple places in your parameterization, HoloPy assumes you want those values tied together (ie they all vary the same and are actually one parameter in the model). So the last Uncertain value around 5 is actually the value for all of the other coordinates.

I think what you probably want to do is give a separate copy of the prior in each place, so that you allow each of those coordinates to vary independently

You can simply type out the prior.Gaussian in each place (which could be a good thing because then you can put a different dither on each one individually), or turn your pg5 into a function:

def pg5():
  return prior.Gaussian(4.97, .08)

then in your sphere code

...
s2 = Sphere(center=(pg5(), pg5(), pg5()) ...
...

That way you will have a separate parameter object in each point and HoloPy will allow the parameters to vary independently.

More about tying parameters in the old fit docs (the inference uses the same parameter engine, we haven't moved all the old docs over) http://holopy.readthedocs.io/en/3.1.1/users/fit_tutorial.html#tying-parameters

[rbnvrw]

Thank you very much for your quick and detailed answer. That is actually a very nice feature that I accidentally stumbled upon, I think it deserves a mention in the documentation!

I have updated and tested my code, it's working now, so I'll close the issue.

[tdimiduk]

If you want to be extra helpful, you could put in a pull request for that documentation upgrade, otherwise someone will get to it eventually.