<p><strong>编辑:重要提示</strong>
这件事已经困扰我一段时间了。我和亚伯拉罕在这里给出的答案是正确的,因为它们增加了x的可变性。然而:请注意,你不能用这种方式简单地增加不确定性来抵消x值中的错误,这样你就回归到“真x”。这个答案中的方法可以告诉你,如果你的x是真的,那么把错误加到x中会如何影响你的回归。如果你的x是错的,这些答案对你没有帮助。x值的误差是一个很难解决的问题,因为它会导致“衰减”和“变量效应误差”。简短的说法是:x中的无偏随机误差会导致回归估计中的<em>偏差</em>。如果你有这个问题,看看卡罗尔,R.J.,鲁伯特,D.,克雷尼西亚努,C.M.和斯特凡斯基,L.A.,2006。<em>非线性模型中的测量误差:现代观点。Chapman和Hall/CRC.,或贝叶斯方法,Gustafson,P.,2003年。<em>统计学和流行病学中的测量误差和错误分类:影响和贝叶斯调整</em>。CRC出版社。最后,我用卡罗尔等人的SIMEX方法和PyMC3解决了我的具体问题。详情见Carstens,H.,Xia,X.和Yadavalli,S.,2017。<em>低成本电能表计量检定校准方法。</em>应用能量,188,第563-575页。在ArXiv上也有</p>
<hr/>
<p>我把上面亚伯拉罕·弗拉克斯曼的答案改成PyMC3,以防有人需要。一些非常小的变化,但可能会混淆无论如何。</p>
<p>第一个是确定性的decorator<code>@Deterministic</code>被一个类似于分布的调用函数<code>var=pymc3.Deterministic()</code>所代替。其次,当生成正态分布随机变量的向量时</p>
<pre><code>rvs = pymc2.rnormal(mu=mu, tau=tau)
</code></pre>
<p>被替换为</p>
<pre><code>rvs = pymc3.Normal('var_name', mu=mu, tau=tau,shape=size(var)).random()
</code></pre>
<p>完整代码如下:</p>
<pre><code>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()
</code></pre>
<p>结果是</p>
<p><a href="https://i.stack.imgur.com/mvhbd.png" rel="nofollow noreferrer">y_error</a></p>
<p>对于x中的错误(请注意变量的“x”后缀):</p>
<pre><code># 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)
</code></pre>
<p>结果是:</p>
<p><a href="https://i.stack.imgur.com/ZOh71.png" rel="nofollow noreferrer">x_error_graph</a></p>
<p>最后的观察是</p>
<pre><code>traceplot(tracex[100:])
plt.tight_layout();
</code></pre>
<p>(结果未显示),我们可以看到<code>sizex</code>似乎由于测量<code>x</code>的误差而遭受“衰减”或“回归稀释”。</p>