与 PyMC3 的贝叶斯相关性

Question

我正在尝试将此 example of Bayesian correlation for PyMC2 转换为 PyMC3，但得到的结果完全不同。最重要的是，多元正态分布的均值很快变为零，而它应该在 400 左右（就像 PyMC2 一样）。因此，估计的相关性很快趋向于 1，这也是错误的。

完整代码可在 notebook for PyMC2 and in this notebook for PyMC3 中找到。

PyMC2的相关代码是

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中如下：

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 版本运行，但显然没有return预期的结果。任何帮助将不胜感激。

Answer 1

pymc.Normal的调用签名是

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。

因此，由于 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 代码生成（类似于）

与 PyMC3 的贝叶斯相关性

Bayesian Correlation with PyMC3

python

correlation

bayesian

pymc

pymc3