python-与PyMC3的贝叶斯相关

我正在尝试将此example of Bayesian correlation for PyMC2转换为PyMC3,但得到的结果完全不同.最重要的是,多元正态分布的均值很快就为零,而它的平均值应为400(与PyMC2一样).因此,估计的相关性迅速变为1,这也是错误的.

完整代码可在此notebook for PyMC2和此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版本可以运行,但显然不会返回预期结果.任何帮助将不胜感激.

解决方法:

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代码就会产生(类似)

标签：bayesian,correlation,pymc,pymc3,python
来源： https://codeday.me/bug/20191026/1937760.html