pymc3: как моделировать коррелированные точки пересечения и наклона в многоуровневой линейной регрессии

В примере Pymc3 для многоуровневой линейной регрессии (пример здесь , с набором данных radon из работы Гельмана и др. (2007)), точки пересечения (для разных округов) и наклоны (для квартир с подвалом и без) имеют нормальный априор. Как я могу смоделировать их вместе с многомерным нормальным априорным, чтобы я мог исследовать корреляцию между ними?

Иерархическая модель, приведенная в примере, выглядит следующим образом:

with pm.Model() as hierarchical_model:
    # Hyperpriors for group nodes
    mu_a = pm.Normal('mu_a', mu=0., sd=100**2)
    sigma_a = pm.HalfCauchy('sigma_a', 5)
    mu_b = pm.Normal('mu_b', mu=0., sd=100**2)
    sigma_b = pm.HalfCauchy('sigma_b', 5)

    # Intercept for each county, distributed around group mean mu_a
    # Above we just set mu and sd to a fixed value while here we
    # plug in a common group distribution for all a and b (which are
    # vectors of length n_counties).
    a = pm.Normal('a', mu=mu_a, sd=sigma_a, shape=n_counties)
    # Intercept for each county, distributed around group mean mu_a
    b = pm.Normal('b', mu=mu_b, sd=sigma_b, shape=n_counties)

    # Model error
    eps = pm.HalfCauchy('eps', 5)

    radon_est = a[county_idx] + b[county_idx] * data.floor.values

    # Data likelihood
    radon_like = pm.Normal('radon_like', mu=radon_est, sd=eps, observed=data.log_radon)
    hierarchical_trace = pm.sample(2000)

И я пытаюсь внести некоторые изменения в приоры

    with pm.Model() as correlation_model:
        # Hyperpriors for group nodes
        mu_a = pm.Normal('mu_a', mu=0., sd=100**2)
        mu_b = pm.Normal('mu_b', mu=0., sd=100**2)

        # here I want to model a and b together
        # I borrowed some code from a multivariate normal model
        # but the code does not work
        sigma = pm.HalfCauchy('sigma', 5, shape=2)

        C_triu = pm.LKJCorr('C_triu', n=2, p=2)
        C = T.fill_diagonal(C_triu[np.zeros((2,2), 'int')], 1)
        cov = pm.Deterministic('cov', T.nlinalg.matrix_dot(sigma, C, sigma))
        tau = pm.Deterministic('tau', T.nlinalg.matrix_inverse(cov))

        a, b = pm.MvNormal('mu', mu=(mu_a, mu_b), tau=tau,
                           shape=(n_counties, n_counties))

        # Model error
        eps = pm.HalfCauchy('eps', 5)
        radon_est = a[county_idx] + b[county_idx] * data.floor.values

        # Data likelihood
        radon_like = pm.Normal('radon_like', mu=radon_est, sd=eps, observed=data.log_radon)
        correlation_trace = pm.sample(2000)

Вот сообщение об ошибке, которое я получил:

  File "<ipython-input-108-ce400c54cc39>", line 14, in <module>
    tau = pm.Deterministic('tau', T.nlinalg.matrix_inverse(cov))

  File "/home/olivier/anaconda3/lib/python3.5/site-packages/theano/gof/op.py", line 611, in __call__
    node = self.make_node(*inputs, **kwargs)

  File "/home/olivier/anaconda3/lib/python3.5/site-packages/theano/tensor/nlinalg.py", line 73, in make_node
    assert x.ndim == 2

AssertionError

Очевидно, я сделал несколько ошибок в матрице ковариации, но я новичок в pymc3 и совершенно новичок в theano, поэтому понятия не имею, как это исправить. Я так понимаю, это должно быть довольно распространенным вариантом использования, так что, может быть, на нем были какие-то примеры? Я просто не могу их найти.

Полный воспроизводимый код и данные можно увидеть на странице примера (ссылка приведена выше). Я не включил его сюда, потому что он слишком длинный, а также я подумал, что те, кто знаком с pymc3, скорее всего, уже хорошо с ним знакомы :)


python pymc3 pymc
person Olivier Ma    schedule 07.09.2016    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Вы забыли добавить одну строку при создании ковариационной матрицы, которую вы пропустили, указав форму MvNormal. Ваша модель должна выглядеть примерно так:

with pm.Model() as correlation_model:
    mu = pm.Normal('mu', mu=0., sd=10, shape=2)
    sigma = pm.HalfCauchy('sigma', 5, shape=2)

    C_triu = pm.LKJCorr('C_triu', n=2, p=2)
    C = tt.fill_diagonal(C_triu[np.zeros((2,2), 'int')], 1.)
    sigma_diag = tt.nlinalg.diag(sigma) # this line
    cov = tt.nlinalg.matrix_dot(sigma_diag, C, sigma_diag)
    tau = tt.nlinalg.matrix_inverse(cov)

    ab = pm.MvNormal('ab', mu=mu, tau=tau, shape=(n_counties, 2))

    eps = pm.HalfCauchy('eps', 5)
    radon_est = ab[:,0][county_idx] + ab[:,1][county_idx] * data.floor.values

    radon_like = pm.Normal('radon_like', mu=radon_est, sd=eps, observed=data.log_radon)
    trace = pm.sample(2000)

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

РЕДАКТИРОВАТЬ

Если вы хотите оценить корреляцию двух переменных апостериорно, вы можете сделать что-то вроде этого.

chain = hierarchical_trace[100:]
x_0 = chain['mu_a']
x_1 = chain['mu_b']
X = np.vstack((x_0, x_1)).T

а затем вы можете запустить следующую модель:

with pm.Model() as correlation:
    mu = pm.Normal('mu', mu=0., sd=10, shape=2)
    sigma = pm.HalfCauchy('sigma', 5, shape=2)

    C_triu = pm.LKJCorr('C_triu', n=2, p=2)
    C = tt.fill_diagonal(C_triu[np.zeros((2,2), 'int')], 1.)
    sigma_diag = tt.nlinalg.diag(sigma)
    cov = tt.nlinalg.matrix_dot(sigma_diag, C, sigma_diag)
    tau = tt.nlinalg.matrix_inverse(cov)

    yl = pm.MvNormal('yl', mu=mu, tau=tau, shape=(2, 2), observed=X)
    trace = pm.sample(5000, pm.Metropolis())

Вы можете заменить x_0 и x_1 в соответствии с вашими потребностями. Например, вы можете захотеть сделать:

x_0 = np.random.normal(chain['mu_a'], chain['sigma_a'])
x_1 = np.random.normal(chain['mu_b'], chain['sigma_b'])
person aloctavodia    schedule 08.09.2016
comment
Спасибо, код работает. Но я уверен, как реализовать ваш альтернативный подход. В hierarchical_model у меня есть среднее значение и стандартное отклонение двух распределений Гаусса, но как они связаны друг с другом? Я предполагаю, что они связаны с наблюдаемыми данными, как вы предположили; но я не уверен, как именно это разыгрывается. correlation_model вычислительно громоздка, но иерархия модели мне ясна. Любые подробности приветствуются! - person Olivier Ma; 09.09.2016