Байесовская корреляция с PyMC3

Question

Я пытаюсь преобразовать этот пример байесовской корреляции для PyMC2 в PyMC3, но получаю совершенно другие результаты. Самое главное, что среднее многовариантного нормального распределения быстро стремится к нулю, тогда как оно должно составлять около 400 (как для PyMC2). Следовательно, расчетная корреляция быстро идет к 1, что также неверно.

Полный код доступен в этом блокноте для PyMC2 и в этом блокноте для 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 работает, но явно не возвращает ожидаемый результат. Любая помощь будет высоко оценен.

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 производит (что-то вроде)