Надежная байесовская корреляция с PyMC3

Question

Это дополнительный вопрос по анализу байесовской корреляции, как описано в этом примере для PyMC2.

Я успешно перенес ненасильный подход, в котором используется многомерное нормальное распределение, на PyMC3, но я борюсь с надежной версией, где вместо этого используется многомерное распределение Student-t.

Проблема в том, что sigma вектор, обозначающий дисперсию, в конечном итоге становится слишком большим, и обратная матрица ковариации больше не может быть вычислена. Следовательно, задняя часть корреляции r почти охватывает весь интервал [-1; 1].

Полный код доступен в этом блокноте для PyMC2 и в этом блокноте для PyMC3.

Соответствующий код для PyMC2:

import pymc
import numpy as np
from scipy.special import gammaln

def analyze(data):
    mu = pymc.Normal('mu', 0, 0.000001, size=2)
    sigma = pymc.Uniform('sigma', 0, 1000, size=2)
    rho = pymc.Uniform('r', -1, 1)

    nu = pymc.Exponential('nu',1/29., 1)
    # we model nu as an Exponential plus one
    @pymc.deterministic
    def nuplus(nu=nu):
        return nu + 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]]))

    # log-likelihood of multivariate t-distribution
    @pymc.stochastic(observed=True)
    def mult_t(value=data.T, mu=mu, tau=precision, nu=nuplus):
        k = float(tau.shape[0])
        res = 0
        for r in value:
            delta = r - mu
            enum1 = gammaln((nu+k)/2.) + 0.5 * np.log(np.linalg.det(tau))
            denom = (k/2.)*np.log(nu*np.pi) + gammaln(nu/2.)
            enum2 = (-(nu+k)/2.) * np.log (1 + (1/nu)*delta.dot(tau).dot(delta.T))
            result = enum1 + enum2 - denom
            res += result[0]
        return res[0,0]

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

Для PyMC3 я могу использовать встроенный MvStudentT распределение, которое ожидает ковариационную матрицу вместо прецизионной матрицы.

import pymc3 as pm
import numpy as np

def mad(data, axis=None):
    return np.median(np.absolute(data - np.median(data, axis)), axis)


def covariance(sigma, rho):
    C = T.alloc(rho, 2, 2)
    C = T.fill_diagonal(C, 1.)
    S = T.diag(sigma)
    return S.dot(C).dot(S)


def analyze_robust(data):
    with pm.Model() as model:
        # priors might be adapted here to be less flat
        mu = pm.Normal('mu', mu=0., tau=0.000001, shape=2, testval=np.median(data, axis=1))
        sigma = pm.Uniform('sigma', lower=0, upper=1000, shape=2, testval=mad(data.T, axis=0))
        rho = pm.Uniform('r', lower=-1., upper=1., testval=0.5)

        cov = pm.Deterministic('cov', covariance(sigma, rho))
        nu = pm.Exponential('nu_minus_one', lam=1./29.) + 1
        mult_t = pm.MvStudentT('mult_t', nu=nu, mu=mu, Sigma=cov, observed=data.T)

    return model

Любые советы, которые могут объяснить поведение под PyMC3, высоко ценятся.

Answer 1

Я обнаружил, что следующее, кажется, дает разумные результаты. Я основал приоры на этой статье, а также использовал ее для проверки корреляции. Вся работа в этой записной книжке.

def mad(data, axis=None):
    return np.median(np.absolute(data - np.median(data, axis)), axis)

def covariance(sigma, rho):
    C = T.alloc(rho, 2, 2)
    C = T.fill_diagonal(C, 1.)
    S = T.diag(sigma)
    return S.dot(C).dot(S)

def analyze_robust(data):
    with pm.Model() as model:
        # priors might be adapted here to be less flat
        mu = pm.Normal('mu', mu=0., sd=100., shape=2, testval=np.median(data.T, axis=1))
        bound_sigma = pm.Bound(pm.Normal, lower=0.)
        sigma = bound_sigma('sigma', mu=0., sd=100., shape=2, testval=mad(data, axis=0))
        rho = pm.Uniform('r', lower=-1., upper=1., testval=0)
        cov = pm.Deterministic('cov', covariance(sigma, rho))
        bound_nu = pm.Bound(pm.Gamma, lower=1.)
        nu = bound_nu('nu', alpha=2, beta=10)
        mult_t = pm.MvStudentT('mult_t', nu=nu, mu=mu, Sigma=cov, observed=data)
    return model