0

在MCMC使用pymc获得最佳拟合值后，我试图制作差异图以测试拟合优度。我的代码如下：

import pymc
import numpy as np
import matplotlib.pyplot as plt, seaborn as sns

# Seeding 
np.random.seed(55555)

# x-data
x = np.linspace(1., 50., 50)

# Gaussian function
def gaus(x, A, x0, sigma): 
        return A*np.exp(-(x-x0)**2/(2*sigma**2))

# y-data
f_true = gaus(x, 10., 25., 10.)
noise = np.random.normal(size=len(f_true)) * 0.2
f = f_true + noise

# y_error
f_err = f*0.05

# Defining the model
def model(x, f):
    A = pymc.Uniform('A', 0., 50., value = 12)
    x0 = pymc.Uniform('x0', 0., 50., value = 20)
    sigma = pymc.Uniform('sigma', 0., 30., value=8)

    @pymc.deterministic(plot=False)
    def gaus(x=x, A=A, x0=x0, sigma=sigma): 
        return A*np.exp(-(x-x0)**2/(2*sigma**2))
    y = pymc.Normal('y', mu=gaus, tau=1.0/f_err**2, value=f, observed=True)
    return locals()

MDL = pymc.MCMC(model(x,f))
MDL.sample(20000, 10000, 1)


# Extract best-fit parameters

A_bf, A_unc = MDL.stats()['A']['mean'], MDL.stats()['A']['standard deviation']
x0_bf, x0_unc = MDL.stats()['x0']['mean'], MDL.stats()['x0']['standard deviation'] 
sigma_bf, sigma_unc = MDL.stats()['sigma']['mean'], MDL.stats()['sigma']['standard deviation']

# Extract and plot results
y_fit = MDL.stats()['gaus']['mean']

plt.clf()
plt.errorbar(x, f, yerr=f_err, color='r', marker='.', label='Observed')
plt.plot(x, y_fit, 'k', ls='-', label='Fit')
plt.legend()
plt.show()

到目前为止一切顺利，并给出了以下情节：使用 MCMC 的最佳拟合图

现在我想使用https://pymc-devs.github.io/pymc/modelchecking.html中第 7.3 节中描述的方法测试拟合优度。为此，我必须先找到 f_sim，所以我在上面几行之后编写了以下代码：

# GOF plot
f_sim = pymc.Normal('f_sim', mu=gaus(x, A_bf, x0_bf, sigma_bf), tau=1.0/f_err**2, size=len(f))
pymc.Matplot.gof_plot(f_sim, f, name='f')
plt.show()

这会给出错误消息AttributeError: 'Normal' object has no attribute 'trace'。我试图在做差异图之前使用 gof_plot 。由于函数的高斯性质，我不认为使用其他分布代替 Normal 是一个好主意。如果有人能让我知道我做错了什么，我将不胜感激。pymc 中的正态分布也没有 Normal_expval 来获得预期值。还有其他方法可以计算 f_exp 吗？谢谢。

4

1 回答 1

0

我意识到 f_sim 实际上是在主要拟合期间定义的 y 值，因为模拟值是蒙特卡罗方法的支柱。所以我提取了最后 10000 次迭代的 y 值并使用 gof_plot 如下：

f_sim = MDL.trace('gaus', chain = None)[:]
pymc.Matplot.gof_plot(f_sim, f, name='f')
plt.show()

现在效果很好！仍然不确定如何获得 f_exp。

于 2015-09-01T19:52:32.690 回答