我试图用高斯(和更复杂)的功能来适应一些数据.我在下面创建了一个小例子.

我的第一个问题是,我做对吗？

我的第二个问题是,如何在x方向添加错误,即在观察/数据的x位置？

很难找到关于如何在pyMC中进行这种回归的好的指南.也许因为它更容易使用一些最小二乘法或类似的方法,但是我最终有很多参数,需要看到我们如何限制它们并比较不同的模型,pyMC似乎是不错的选择.

import pymc

import numpy as np

import matplotlib.pyplot as plt; plt.ion()

x = np.arange(5,400,10)*1e3

# Parameters for gaussian

amp_true = 0.2

size_true = 1.8

ps_true = 0.1

# Gaussian function

gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps

f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )

# add noise to the data points

noise = np.random.normal(size=len(x)) * .02

f = f_true + noise

f_error = np.ones_like(f_true)*0.05*f.max()

# define the model/function to be fitted.

def model(x, f):

amp = pymc.Uniform('amp', 0.05, 0.4, value= 0.15)

size = pymc.Uniform('size', 0.5, 2.5, value= 1.0)

ps = pymc.Normal('ps', 0.13, 40, value=0.15)

@pymc.deterministic(plot=False)

def gauss(x=x, amp=amp, size=size, ps=ps):

e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.))

return amp*np.exp(e)+ps

y = pymc.Normal('y', mu=gauss, tau=1.0/f_error**2, value=f, observed=True)

return locals()

MDL = pymc.MCMC(model(x,f))

MDL.sample(1e4)

# extract and plot results

y_min = MDL.stats()['gauss']['quantiles'][2.5]

y_max = MDL.stats()['gauss']['quantiles'][97.5]

y_fit = MDL.stats()['gauss']['mean']

plt.plot(x,f_true,'b', marker='None', ls='-', lw=1, label='True')

plt.errorbar(x,f,yerr=f_error, color='r', marker='.', ls='None', label='Observed')

plt.plot(x,y_fit,'k', marker='+', ls='None', ms=5, mew=2, label='Fit')

plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5)

plt.legend()

我意识到,我可能需要运行更多的迭代,最后使用刻录和稀释.绘制数据和拟合的图在下面看到.

pymc.Matplot.plot(MDL)数字看起来像这样,显示出很好的峰值分布.这很好,对吧？