Байесовские методы гетероскедастичности

Это будет короткий пост, потому что знаете что? Вот только не надо, чтобы он был длинным. Иногда нужно просто что-то увидеть, и ты сразу это понимаешь.

Иногда ваши данные гетероскедастичны, и вы не хотите делать глупости, чтобы с этим справиться. Байесовские методы спасают в обоих случаях.

Гетероскедастичность — это, в основном, когда дисперсия вашей целевой переменной отличается для заданных значений ваших данных. Или, выражаясь более математически, дисперсия y зависит от x.

Сначала давайте создадим несколько поддельных данных, чтобы вы знали, с чем мы имеем дело.

Код стандартной линейной модели в PyMC3 приведен ниже, и он будет эталоном.

with pm.Model() as simple_model:
    # Priors
    m = pm.Normal('m', mu=0, sd=100**2)
    b = pm.Normal('b', mu=0, sd=100**2)
    sig = pm.HalfNormal('sigma', sd=100**2)
    
    # Likelihood    
    rsq = pm.Deterministic('rsq', r2_pymc3(y, m*x + b))
    
    obs = pm.Normal('y', mu=m*x + b, sd=sig, observed=y)
    
    trace_simple_model = pm.sample(2000)
    
    pm.traceplot(trace_simple_model)
    graph.show()

На самом деле имеет дело с гетероскедастичностью

Метод 1. Коши как функция правдоподобия

# Method 1. Cauchy Method
with pm.Model() as cauchy_model:
    # Priors
    m = pm.Normal('m', mu=0, sd=100**2)
    b = pm.Normal('b', mu=0, sd=100**2)
    sig = pm.HalfNormal('sigma', sd=100**2)
    
    # Likelihood    
    rsq = pm.Deterministic('rsq', r2_pymc3(y, m*x + b))
    
    obs = pm.Cauchy('y', alpha=m*x + b, beta=sig, observed=y)
    
    trace_cauchy_model = pm.sample(2000)
    
    pm.traceplot(trace_cauchy_model)
    graph.show()

Обратите внимание, что дисперсия m и b более узкая, чем в обычной байесовской регрессии. Однако точка пересечения слишком высока, а наклон слишком низок, как вы можете видеть на апостериорных графиках ниже.

Метод 2. Моделирование дисперсии по данным

# Method 2. modeling 
with pm.Model() as fancy_model:
    n_samples = len(y)
    data, target = x, y
    
    # Priors
    m = pm.Normal('m', mu=0, sd=100**2)
    b = pm.Normal('b', mu=0, sd=100**2)
    
    sig_m = pm.Normal('sigma_m', mu=0, sd=100**2)
    sig_b = pm.Normal('sigma_b', mu=0, sd=100**2)
    sig_given_data = pm.Deterministic('sig', pm.math.log(1 + pm.math.exp(sig_m * data + sig_b)))
    
    # Likelihood    
    rsq = pm.Deterministic('rsq', r2_pymc3(y, m*x + b))
    
    obs = pm.Normal('y', mu=m*data + b, sd=sig_given_data, observed=target)
    
    trace_fancy_model = pm.sample(2000)
    
    pm.traceplot(trace_fancy_model)
    graph.show()

Обратите внимание, что дисперсия m и b все еще более узкая, чем в обычной байесовской регрессии. Однако, в отличие от распределения Коши, наклон и точка пересечения оцениваются лучше.

Вывод

Байесовские методы прекрасны. Так что используйте их.

Полный и более подробный блокнот Jupyter можно найти здесь.