使用pyMCMC/pyMC拟合非线性函数到数据/观测值
我正在尝试用高斯函数（还有更复杂的函数）来拟合一些数据。我在下面创建了一个小例子。
我第一个问题是，我这样做对吗？
我第二个问题是，我怎么在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)生成的图看起来像这样，显示了很好分布的峰。这是好的，对吧？
2 个回答
编辑：重要说明
这件事让我困扰了一段时间。我和Abraham在这里给出的答案是正确的，因为它们为x增加了变化性。然而，请注意，你不能仅仅通过这种方式添加不确定性来抵消x值中的错误，以便对“真实的x”进行回归。这个答案中的方法可以告诉你如果你有真实的x，添加错误到x会如何影响你的回归。如果你的x测量不准确，这些答案就帮不了你。x值中的错误是一个非常棘手的问题，因为它会导致“衰减”和“变量中的错误效应”。简单来说，x中存在无偏随机错误会导致你的回归估计出现偏差。如果你遇到这个问题，可以参考Carroll, R.J., Ruppert, D., Crainiceanu, C.M. 和 Stefanski, L.A.在2006年出版的非线性模型中的测量误差：现代视角。Chapman和Hall/CRC出版，或者如果你想用贝叶斯方法，可以看Gustafson, P.在2003年出版的统计学和流行病学中的测量误差和误分类：影响和贝叶斯调整。CRC出版社。我最终使用Carroll等人的SIMEX方法和PyMC3解决了我的具体问题。详细信息见Carstens, H., Xia, X. 和 Yadavalli, S.在2017年发表的低成本能量计校准方法用于测量和验证。《应用能源》，188，563-575页。该文也可以在ArXiv上找到。
我把Abraham Flaxman上面的答案转换成了PyMC3，以防有人需要。做了一些非常小的修改，但仍然可能会让人困惑。
首先，确定性装饰器@Deterministic被一个类似分布的调用函数var=pymc3.Deterministic()替代。其次，在生成一个正态分布随机变量的向量时，
rvs = pymc2.rnormal(mu=mu, tau=tau)
被替换为
rvs = pymc3.Normal('var_name', mu=mu, tau=tau,shape=size(var)).random()
完整代码如下：
import numpy as np
from pymc3 import *
import matplotlib.pyplot as plt
# set random seed for reproducibility
np.random.seed(12345)
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()
with Model() as model3:
amp = Uniform('amp', 0.05, 0.4, testval= 0.15)
size = Uniform('size', 0.5, 2.5, testval= 1.0)
ps = Normal('ps', 0.13, 40, testval=0.15)
gauss=Deterministic('gauss',amp*np.exp(-1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.)))+ps)
y =Normal('y', mu=gauss, tau=1.0/f_error**2, observed=f)
start=find_MAP()
step=NUTS()
trace=sample(2000,start=start)
# extract and plot results
y_min = np.percentile(trace.gauss,2.5,axis=0)
y_max = np.percentile(trace.gauss,97.5,axis=0)
y_fit = np.percentile(trace.gauss,50,axis=0)
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=1, label='Fit')
plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5)
plt.legend()
这将产生
对于x中的错误（注意变量的'x'后缀）：
# define the model/function to be fitted in PyMC3:
with Model() as modelx:
x_obsx = pm3.Normal('x_obsx',mu=x, tau=(1e4)**-2, shape=40)
ampx = Uniform('ampx', 0.05, 0.4, testval=0.15)
sizex = Uniform('sizex', 0.5, 2.5, testval=1.0)
psx = Normal('psx', 0.13, 40, testval=0.15)
x_pred = Normal('x_pred', mu=x_obsx, tau=(1e4)**-2*np.ones_like(x_obsx),testval=5*np.ones_like(x_obsx),shape=40) # this allows error in x_obs
gauss=Deterministic('gauss',ampx*np.exp(-1*(np.pi**2*sizex*x_pred/(3600.*180.))**2/(4.*np.log(2.)))+psx)
y = Normal('y', mu=gauss, tau=1.0/f_error**2, observed=f)
start=find_MAP()
step=NUTS()
tracex=sample(20000,start=start)
这将产生：
最后一个观察是，当进行
traceplot(tracex[100:])
plt.tight_layout();
（结果未显示），我们可以看到sizex似乎因为x的测量误差而遭受了“衰减”或“回归稀释”。
我第一个问题是，我这样做对吗？
对的！你需要加入一个“烧入期”，你已经知道了。我通常会把前一半的数据样本丢掉。你不需要做任何的“稀疏处理”，但有时候这样可以让你后面的工作处理得更快，存储得更小。
我唯一建议的就是设置一个随机种子，这样你的结果才能“可重复”：np.random.seed(12345) 这样就可以了。
哦，如果我真的要多说一点，我会建议你使用import seaborn，这样可以让matplotlib的结果看起来更漂亮。
我第二个问题是，如何在x方向上添加误差，也就是在观察数据的x位置上？
一种方法是为每个误差添加一个潜在变量。在你的例子中这样做是可行的，但如果你有很多观察数据，这样做就不太现实了。我给你一个小例子，帮助你入门：
# add noise to observed x values
x_obs = pm.rnormal(mu=x, tau=(1e4)**-2)
# define the model/function to be fitted.
def model(x_obs, f):
amp = pm.Uniform('amp', 0.05, 0.4, value= 0.15)
size = pm.Uniform('size', 0.5, 2.5, value= 1.0)
ps = pm.Normal('ps', 0.13, 40, value=0.15)
x_pred = pm.Normal('x', mu=x_obs, tau=(1e4)**-2) # this allows error in x_obs
@pm.deterministic(plot=False)
def gauss(x=x_pred, 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 = pm.Normal('y', mu=gauss, tau=1.0/f_error**2, value=f, observed=True)
return locals()
MDL = pm.MCMC(model(x_obs, f))
MDL.use_step_method(pm.AdaptiveMetropolis, MDL.x_pred) # use AdaptiveMetropolis to "learn" how to step
MDL.sample(200000, 100000, 10) # run chain longer since there are more dimensions
看起来如果你的x和y都有噪声，得到好的答案可能会很困难：
这里有一个收集这些内容的笔记本。