0

我正在尝试使用 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.
4

2 回答 2

0

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

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

我知道您可以使用一种解决方法，即检查mu变量是否为 a 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 问题跟踪器中提交错误报告。

于 2014-08-20T05:47:16.210 回答
0

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

flux = mc.Normal('flux',mu=gaussian.value,tau=err_map,value=arr,observed=True,doc='Observed Flux')
于 2016-03-07T19:19:57.997 回答