Как рассчитать симулированные значения при построении графика несоответствия для хорошего соответствия?

Question

Я пытаюсь сделать график несоответствия для тестирования доброкачественности после получения наилучших значений соответствия 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() 

До сих пор так хорошо и дает следующий сюжет:

Теперь я хочу, чтобы проверить благость-из-приступе с использованием метода, описанного в разделе 7.3 в https://pymc-devs.github.io/pymc/modelchecking.html. Для этого я должен найти 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» не имеет атрибута «след». Я пытаюсь использовать gof_plot перед выполнением графика несоответствия. Я не думаю, что использование другого распределения вместо Normal было бы хорошей идеей из-за гауссовского характера функции. Я был бы очень признателен, если бы кто-нибудь мог сообщить мне, что я делаю неправильно. Также нормальное распределение в pymc не имеет Normal_expval для получения ожидаемых значений. Есть ли другой способ вычисления f_exp? Благодарю.

Answer 1

Я понял, что f_sim - это фактически значения y, определенные во время основной подгонки, поскольку имитируемые значения являются основой метода montecarlo. Таким образом, я извлек Y значения за последние 10000 итераций и используется gof_plot следующим образом:

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

теперь работает отлично! Все еще не уверен, как получить f_exp.