2016-10-16 11 views
3

私はこのexample of Bayesian correlation for PyMC2をPyMC3に変換しようとしていますが、全く異なる結果になります。最も重要な点は、多変量正規分布の平均は急速にゼロになりますが、PyMC2の場合と同様に約400になるはずです。その結果、推定された相関は速やかに1に向かい、これも間違っている。Bayesian PyMC3との相関

完全なコードはnotebook for PyMC2notebook for PyMC3にあります。次のように

PyMC2に関連するコードがPyMC3に上記のコードの

def analyze(data): 
    # priors might be adapted here to be less flat 
    mu = pymc.Normal('mu', 0, 0.000001, size=2) 
    sigma = pymc.Uniform('sigma', 0, 1000, size=2) 
    rho = pymc.Uniform('r', -1, 1) 

    @pymc.deterministic 
    def precision(sigma=sigma,rho=rho): 
     ss1 = float(sigma[0] * sigma[0]) 
     ss2 = float(sigma[1] * sigma[1]) 
     rss = float(rho * sigma[0] * sigma[1]) 
     return np.linalg.inv(np.mat([[ss1, rss], [rss, ss2]])) 

    mult_n = pymc.MvNormal('mult_n', mu=mu, tau=precision, value=data.T, observed=True) 

    model = pymc.MCMC(locals()) 
    model.sample(50000,25000) 

マイポートがあるさ:

def precision(sigma, rho): 
    C = T.alloc(rho, 2, 2) 
    C = T.fill_diagonal(C, 1.) 
    S = T.diag(sigma) 
    return T.nlinalg.matrix_inverse(T.nlinalg.matrix_dot(S, C, S)) 


def analyze(data): 
    with pm.Model() as model: 
     # priors might be adapted here to be less flat 
     mu = pm.Normal('mu', mu=0., sd=0.000001, shape=2, testval=np.mean(data, axis=1)) 
     sigma = pm.Uniform('sigma', lower=1e-6, upper=1000., shape=2, testval=np.std(data, axis=1)) 
     rho = pm.Uniform('r', lower=-1., upper=1., testval=0) 

     prec = pm.Deterministic('prec', precision(sigma, rho)) 
     mult_n = pm.MvNormal('mult_n', mu=mu, tau=prec, observed=data.T) 

    return model 

model = analyze(data) 
with model: 
    trace = pm.sample(50000, tune=25000, step=pm.Metropolis()) 

PyMC3のバージョンが実行されますが、はっきりと期待される結果を返しません。どんな助けも高く評価されます。

答えて

3

pymc.Normalのコールサインはpymc.Normalの第三の位置の引数は、標準偏差、sdtauないこと

In [125]: pymc.Normal? 
Init signature: pymc.Normal(self, *args, **kwds) 
Docstring: 
N = Normal(name, mu, tau, value=None, observed=False, size=1, trace=True, rseed=True, doc=None, verbose=-1, debug=False) 

注意あります。したがって

、対応pymc3コードはtau = 1/sigma**2ため

mu = pm.Normal('mu', mu=0., tau=0.000001, shape=2, ...) 

又は

mu = pm.Normal('mu', mu=0., sd=math.sqrt(1/0.000001), shape=2, ...) 

を使用する必要が

mu = Normal('mu', 0, 0.000001, size=2) 

pymcコードを使用するため。この1回の変更で


、あなたのpymc3コードは(のようなもの)を生成

enter image description here

+0

愚かな間違い、予想通りミューための先行がmisspecifiedされている場合、それが動作しないという意味になります。これを指摘してくれてありがとう。 – sebp

関連する問題