Pymc 2D高斯拟合
我正在尝试用pymc将一个预定义的二维高斯函数拟合到一些观察到的数据上。但是我一直遇到错误，最近一次出现的错误是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.
2 个回答
0
@mc.deterministic这个装饰器会返回一个叫做deterministic的变量。要获取这个变量的值，可以使用value这个属性。
flux = mc.Normal('flux',mu=gaussian.value,tau=err_map,value=arr,observed=True,doc='Observed Flux')
0
我之前也遇到过类似的问题，但从来没有机会找到根本原因。你代码中有问题的那一行是关于观察到的 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报告。