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

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

Теперь я хочу проверить соответствие, используя метод, описанный в разделе 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: «Нормальный» объект не имеет атрибута «трассировка». Я пытаюсь использовать gof_plot, прежде чем делать график несоответствия. Я не думаю, что использование другого распределения вместо нормального было бы хорошей идеей из-за гауссовского характера функции. Я был бы очень признателен, если бы кто-нибудь мог сообщить мне, что я делаю неправильно. Также нормальное распределение в pymc не имеет Normal_expval для получения ожидаемых значений. Есть ли другой способ рассчитать f_exp? Спасибо.


python mcmc goodness-of-fit pymc
person Silentrash    schedule 31.08.2015    source источник
comment
На некоторые из ваших вопросов можно ответить по адресу stackoverflow.com/questions/30731681/.   -  person    schedule 01.09.2015
comment
Я смотрел это раньше, но я не был полностью уверен, как вписать это описание в контекст моей проблемы. Я думаю, что, наконец, понял это. Я опубликую свое решение сейчас. Спасибо, в любом случае.   -  person Silentrash    schedule 01.09.2015

Ответы (1)


arrow_upward
0
arrow_downward

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

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

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

person Silentrash    schedule 01.09.2015