0

pymc を使用して MCMC で最適な適合値を取得した後、適合度をテストするための不一致プロットを作成しようとしています。私のコードは次のようになります：

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 を使用しようとしています。関数のガウス特性のため、通常の代わりに他の分布を使用することは良い考えではないと思います。誰かが私が間違っていることを教えてくれれば、本当に感謝しています。また、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 に答える