与PyMC3的贝叶斯相关 [英] Bayesian Correlation with PyMC3

问题描述

我正在尝试将此将PyMC2的贝叶斯相关性示例转换为PyMC3，但是得到完全不同的结果.最重要的是，多元正态分布的均值迅速变为零，而它的平均值应为400(与PyMC2一样).因此，估计的相关性迅速变为1，这也是错误的.

I'm trying to convert this example of Bayesian correlation for PyMC2 to PyMC3, but get completely different results. Most importantly, the mean of the multivariate Normal distribution quickly goes to zero, whereas it should be around 400 (as it is for PyMC2). Consequently, the estimated correlation quickly goes towards 1, which is wrong as well.

完整代码可在此 PyMC2笔记本中获得，并在此 PyMC3笔记本中.

The full code is available in this notebook for PyMC2 and in this notebook for PyMC3.

PyMC2的相关代码是

The relevant code for PyMC2 is

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) 

我上面的代码到PyMC3的端口如下:

My port of the above code to PyMC3 is as follows:

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版本可以运行，但是显然不会返回预期的结果.任何帮助将不胜感激.

The PyMC3 version runs, but clearly does not return the expected result. Any help would be highly appreciated.

推荐答案

pymc.Normal的呼叫签名为

The call signature of pymc.Normal is

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)

请注意，pymc.Normal的第三个位置参数是tau，而不是标准偏差sd.

Notice that the third positional argument of pymc.Normal is tau, not the standard deviation, sd.

因此，由于pymc代码使用

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

应使用相应的pymc3代码

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, ...)

自tau = 1/sigma**2起.

通过这一更改，您的pymc3代码会生成(类似)

With this one change, your pymc3 code produces (something like)

这篇关于与PyMC3的贝叶斯相关的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！