Pymc 2d gaussian fitting

Question

I am trying to fit a predefined 2d gaussian function to some observed data with pymc. I keep running into errors and the last one I got was ValueError: setting an array element with a sequence. I understand what the error means, but I am not sure where the error is occurring in the code. My naive guess would be the random variables are being set to some array elements. Any suggestions would be much appreciated. Here is my code so far:

import pymc as mc
import numpy as np
import pyfits as pf

arr = pf.getdata('img.fits')
x=y=np.arange(0,71)
xx,yy=np.meshgrid(x,y)
err_map = pf.getdata('imgwht.fits')
def model((x,y),arr):
    amp = mc.Uniform('amp',lower=-1,upper=1,doc='Amplitude')
    x0 = mc.Uniform('x0',lower=21,upper=51,doc='xo')
    y0 = mc.Uniform('y0',lower=21,upper=51,doc='yo')
    sigx = mc.Uniform('sigx',lower=0.1,upper=10,doc='Sigma in X')
    sigy = mc.Uniform('sigy',lower=0.1,upper=10,doc='Sigma in Y')
    thta = mc.Uniform('theta',lower=0,upper=2*np.pi,doc='Rotation')
    os = mc.Uniform('c',lower=-1,upper=1,doc='Vertical offset')

    @mc.deterministic(plot=False,trace=False)
    def gaussian((x, y)=(xx,yy), amplitude=amp, xo=x0, yo=y0, sigma_x=sigx, sigma_y=sigy, theta=thta, offset=os):
        xo = float(xo)
        yo = float(yo)    
        a = (mc.cos(theta)**2)/(2*sigma_x**2) + (mc.sin(theta)**2)/(2*sigma_y**2)
        b = -(mc.sin(2*theta))/(4*sigma_x**2) + (mc.sin(2*theta))/(4*sigma_y**2)
        c = (mc.sin(theta)**2)/(2*sigma_x**2) + (mc.cos(theta)**2)/(2*sigma_y**2)
        gauss = offset+amplitude*mc.exp(-1*(a*((x-xo)**2)+2*b*(x-xo)*(y-yo)+c*((y-yo)**2)))
        return gauss
    flux = mc.Normal('flux',mu=gaussian,tau=err_map,value=arr,observed=True,doc='Observed Flux')
    return locals()
mdl = mc.MCMC(model((xx,yy),arr))
mdl.sample(iter=1e5,burn=9e4)

Full traceback:

File "model.py", line 31, in <module>
    mdl = mc.MCMC(model((xx,yy),arr))
File "model.py", line 29, in model
    flux = mc.Normal('flux',mu=gaussian,tau=err_map,value=arr,observed=True,doc='Observed Flux')
File "/usr/lib64/python2.7/site-packages/pymc/distributions.py", line 318, in __init__
**arg_dict_out)
File "/usr/lib64/python2.7/site-packages/pymc/PyMCObjects.py", line 761, in __init__
verbose=verbose)
File "/usr/lib64/python2.7/site-packages/pymc/Node.py", line 219, in __init__
Node.__init__(self, doc, name, parents, cache_depth, verbose=verbose)
File "/usr/lib64/python2.7/site-packages/pymc/Node.py", line 129, in __init__
self.parents = parents
File "/usr/lib64/python2.7/site-packages/pymc/Node.py", line 152, in _set_parents
self.gen_lazy_function()
File "/usr/lib64/python2.7/site-packages/pymc/PyMCObjects.py", line 810, in gen_lazy_function
self._logp.force_compute()
File "LazyFunction.pyx", line 257, in pymc.LazyFunction.LazyFunction.force_compute (pymc/LazyFunction.c:2409)
File "/usr/lib64/python2.7/site-packages/pymc/distributions.py", line 2977, in wrapper
return f(value, **kwds)
File "/usr/lib64/python2.7/site-packages/pymc/distributions.py", line 2168, in normal_like
return flib.normal(x, mu, tau)
ValueError: setting an array element with a sequence.

Welcome to Stack Overflow, please take the [Tour](http://stackoverflow.com/tour). Questions seeking debugging help ("why isn't this code working?") must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Questions without a clear problem statement are not useful to other readers. See [How to create a Minimal, Complete, and Verifiable example](http://stackoverflow.com/help/mcve). — DavidPostill, Aug 16 '14 at 18:58
I'm getting `IOError: File does not exist: 'img.fits'`; can you link to the `.fits` files that your code requires? — Abraham D Flaxman, Aug 18 '14 at 06:28
@AbrahamDFlaxman Here's the link to img.fits: http://s000.tinyupload.com/?file_id=65237298882568362358 and for imgwht.fits: http://s000.tinyupload.com/?file_id=02197311909809956804 — bmsmith, Aug 18 '14 at 21:23

score 0 · Answer 1 · answered Aug 20 '14 at 05:47

I've run into an issue like this before, but never had a chance to track it down to its source. The problem line in your code is the one for the observed Stochastic:

flux = mc.Normal('flux',mu=gaussian,tau=err_map,value=arr,observed=True,doc='Observed Flux')

I know a work-around that you can use, which is to check if the mu variable is a pymc.Node, and only find the likelihood if it is not:

@mc.observed
def flux(mu=gaussian,tau=err_map,value=arr):
    if isinstance(mu, mc.Node):
        return 0
    else:
        return mc.normal_like(value, mu, tau)

I think it would be worth filing a bug report in the PyMC github issue tracker if you have time.

score 0 · Answer 2 · answered Mar 07 '16 at 19:19

0

The @mc.deterministic decorator returns a deterministic variable. To get the value of the variable, use the attribute value.

flux = mc.Normal('flux',mu=gaussian.value,tau=err_map,value=arr,observed=True,doc='Observed Flux')

answered Mar 07 '16 at 19:19

Lila

509
4
5

Pymc 2d gaussian fitting

2 Answers2