Соответствие logp в PyMC3 - PullRequest
Новая
       2

Соответствие logp в PyMC3

0 голосов
/ 14 февраля 2020

Я новичок в PyMC3 и проверяю свое понимание. У меня есть простая линейная модель y = a + b * x + error. 

from pymc3 import Model, HalfCauchy, Normal, sample
from pymc3 import Lognormal
from scipy.stats import norm, lognorm

import numpy as np

# Parameters
# One observation, linear model: y = a + b*x + e
a_prior_mu, a_prior_sg = 0, 5
b_prior_mu, b_prior_sg = 0.5, 0.1
sg_sg_prior = 3
a, b = 0.1, 0.55
sg_y = 1.5

# Data
size = 1
np.random.seed = 1234
x = np.random.randn(size)
e = np.random.randn(size)
y_obs = a + b * x + e

with Model() as model_wo_Cauchy:
    a = Normal('a', mu=a_prior_mu, sigma=a_prior_sg)
    b = Normal('b', mu=b_prior_mu, sigma=b_prior_sg)
    likelihood = Normal('y', mu=a + b*x, sigma=sg_y, observed=y_obs)

with Model() as model_wi_Cauchy:
    a = Normal('a', mu=a_prior_mu, sigma=a_prior_sg)
    b = Normal('b', mu=b_prior_mu, sigma=b_prior_sg)
    sg = Lognormal('sg',  mu=0, sigma=sg_sg_prior)
    likelihood = Normal('y', mu=a + b*x, sigma=sg, observed=y_obs)

Я пытаюсь объяснить вероятность: 

# Example of posterior
a_hat, b_hat, sg_hat = 2, 0.6, 2

# Explain log-likelihood
dic_tmp_wo_Cauchy = {model_wo_Cauchy.free_RVs[0]: a_hat,
                     model_wo_Cauchy.free_RVs[1]: b_hat}
dic_tmp_wi_Cauchy = {model_wi_Cauchy.free_RVs[0]: a_hat,
                     model_wi_Cauchy.free_RVs[1]: b_hat,
                     model_wi_Cauchy.free_RVs[2]: np.log(sg_hat)}

(Простите за наименование, я сначала попробовал с HalfCauchy, затем переключился на Lognormal, чтобы посмотреть, получу ли я лучшие результаты)

Репликация правдоподобия журналов: 

# Check
log_p_a = norm.logpdf(a_hat, a_prior_mu, a_prior_sg)
log_p_b = norm.logpdf(b_hat, b_prior_mu, b_prior_sg)
log_p_s = lognorm.logpdf(sg_hat, sg_sg_prior)
log_p_yGab = norm.logpdf(y_obs, a_hat + b_hat * x, sg_y)[0]
log_p_yGabs = norm.logpdf(y_obs, a_hat + b_hat * x, sg_hat)[0]

# Comparison
print(model_wo_Cauchy.logp(dic_tmp_wo_Cauchy))    # -3.3623336246196764
print(log_p_yGab + log_p_a + log_p_b)             # -3.3623336246196764
print(model_wi_Cauchy.logp(dic_tmp_wi_Cauchy))    # -5.557233310413178
print(log_p_yGabs + log_p_a + log_p_b + log_p_s)  # -6.250380490973123

Без логнормальной части распределения вещи совпадают, но я не могу сопоставить их с логнормальный Как ни странно, если наблюдение sg_hat было = 1, на самом деле результаты совпадают с последними 2 утверждениями, что приводит к -5.3659194602060944).

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

sg.logp({'sg': sg_hat})

Но я получаю ошибку: 

sg.logp({'sg': sg_hat})
Traceback (most recent call last):

  File "<ipython-input-13-0c061388ccc9>", line 1, in <module>
    sg.logp({'sg': sg_hat})

AttributeError: 'TransformedRV' object has no attribute 'logp'

В терминах ошибки, копаясь вокруг, я нашел: https://datascience.stackexchange.com/questions/62970/use-of-pymc-distribution-dist I, который я использовал для проверки этого: 

with Model() as model_chk:
    sg_chk = Lognormal('sg_chk',  mu=0, sigma=sg_sg_prior)    # -2.0442426559793496
print(model_chk.logp({'sg_chk_log__': np.log(sg_hat)}))       # -2.7373898365392946
print(log_p_s)

Однако, если я делаю это: 

with Model() as model_chk_no_transf:
    sg_chk = Lognormal('sg_chk',  mu=0, sigma=sg_sg_prior, transform=None)
print(model_chk_no_transf.logp({'sg_chk': sg_hat}))           # -2.737389836539295
print(log_p_s)                                                # -2.7373898365392946

Из этого я делаю вывод, что я вычисляю априор в другом чем PyMC3, то есть я бы делал эквивалент не-преобразования, что звучит странно.

Любое понимание было бы очень полезно!

Большое спасибо

...