二维高斯拟合

Question

我试图用pymc将预定义的2d高斯函数拟合到一些观察到的数据。我不断遇到错误，我得到的最后一个是ValueError: setting an array element with a sequence.，我理解错误的意思，但我不确定错误在代码中发生的地方。我天真的猜测是随机变量被设置为一些数组元素。如有任何建议，我将不胜感激。下面是我到目前为止的代码:

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)

完整回溯:

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.

Answer 1

我以前遇到过这样的问题，但从来没有机会追踪到它的来源。代码中的问题行是观察到的Stochastic:

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

我知道一个可以使用的解决方法，这是检查mu变量是否为pymc.Node，并且仅在它不是时才找到可能性:

@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)

我认为如果你有时间的话，在PyMC github问题跟踪器中提交一个bug报告是值得的。

Answer 2

@mc.deterministic装饰器返回一个deterministic变量。要获取变量的值，请使用属性value。

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