0

私はpymc2に次のモデルを持っています:

import pymc
from scipy.stats import gamma

alpha = pymc.Uniform('alpha', 0.01, 2.0)
scale = pymc.Uniform('scale', 1.0, 4.0)

@pymc.deterministic(plot=False)
def beta(scale=scale):
    return 1.0 / scale

@pymc.potential
def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)):
    dist = gamma(alpha, loc=0., scale=scale)
    fp = 1.0 - dist.cdf(lmin)
    return -(n+1)*np.log(fp)

obs = pymc.Gamma("obs", alpha=alpha, beta=beta, value=sample, observed=True)

このモデルの物理的背景は、銀河の光度関数(LF)、つまり銀河が光度 L を持つ確率です。銀河の種類によっては、LF は単なるガンマ関数です。銀河の調査では通常、ターゲットのかなりの部分、特に光度の低いターゲットを見逃しているため、データの切り捨ての可能性が説明されています。このモデルでは、以下のすべてが欠けていますlmin

この方法の詳細は、ケリーらによるこの論文で見つけることができます。

このモデルは機能します。モデルを実行するMAPと、シミュレートされたデータからパラメーターをMCMC回復できますが、成長するにつれて不確実性が高まります。alphascalesamplelmin

ここで、ガウス測定誤差を挿入したいと思います。簡単にするために、すべてのデータの精度は同じです。エラーも含める可能性を変更していません。

alpha = pymc.Uniform('alpha', 0.01, 2.0)
scale = pymc.Uniform('scale',1.0, 4.0)
sig = 0.1
tau = math.pow(sig, -2.0)  

@pymc.deterministic(plot=False)
def beta(scale=scale):
    return 1.0 / scale

@pymc.potential
def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)):
    dist = gamma(alpha, loc=0., scale=scale)
    fp = 1.0 - dist.cdf(lmin)
    return -(n+1) * np.log(fp)

dist = pymc.Gamma("dist", alpha=alpha, beta=beta)
obs = pymc.Normal("obs", mu=dist, tau=tau, value=sample, observed=True)

しかし、このモデルは機能しないため、ここで何か間違ったことをしていることは確かです。このモデルで実行すると、とpymc.MAPの初期値が回復しますalphascale

vals = {'alpha': alpha, 'scale': scale, 'beta': beta, 
   'p_factor': p_factor, 'obs': obs, 'dist': dist}
M2 = pymc.MAP(vals)
M2.fit()
print M2.alpha.value, M2.scale.value
>>> (array(0.010000000006018368), array(1.000000000833973))

を実行するpymc.MCMCとalpha、betaまったくトレースされません。

M = pymc.MCMC(vals)
M.sample(10000, burn=5000)
...
M.stats()['alpha']
>>> {'95% HPD interval': array([ 0.01000001,  0.01000502]),
'mc error': 2.1442678276712383e-07,
'mean': 0.010001588137798096,
'n': 5000,
'quantiles': {2.5: 0.0100000088679046,
25: 0.010000382359859467,
50: 0.010001100377476166,
75: 0.010001668672799679,
97.5: 0.0100050194240779},
'standard deviation': 2.189828287191421e-06}

再び初期値。実際alpha、たとえば 0.02 で開始するように変更すると、復元された値alphaは 0.02 になります。

これは、作業モデルとシミュレートされたデータを含むノートブックです。

これは、エラー モデルとシミュレートされたデータを含むノートブックです。

この作業を行うためのガイダンスは本当にありがたいです。

4

1 に答える 1

1

変更するだけで十分なようです

dist = pymc.Gamma("dist", alpha=alpha, beta=beta)

に

dist = pymc.Gamma("dist", alpha=alpha, beta=beta, value=sample)

サンプリングされたデータは の妥当な初期値ですdist。alphaとにかく、他の初期値（ゼロの配列など）がサンプリングしないという問題を再び引き起こすため、ロジックがわかりませんbeta。

于 2014-03-07T00:06:51.803 に答える