Многомерная линейная регрессия в pymc3

Question

Я недавно начал учиться pymc3 после исключительно использования emcee целую вечность, и я сталкиваюсь с некоторыми концептуальными проблемами.

Я практикуюсь в главе 7 " Подгонка модели к данным Хогга". Это подразумевает подгонку mcmc к прямой линии с произвольными 2-мерными неопределенностями. Я сделал это довольно легко в emcee, но pymc дает мне некоторые проблемы.

По сути, это сводится к использованию многомерной гауссовской вероятности.

Вот что у меня так далеко.

from pymc3 import  *

import numpy as np
import matplotlib.pyplot as plt

size = 200
true_intercept = 1
true_slope = 2

true_x = np.linspace(0, 1, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * true_x
# add noise

# here the errors are all the same but the real world they are usually not!
std_y, std_x = 0.1, 0.1 
y = true_regression_line + np.random.normal(scale=std_y, size=size)
x = true_x + np.random.normal(scale=std_x, size=size)

y_err = np.ones_like(y) * std_y
x_err = np.ones_like(x) * std_x

data = dict(x=x, y=y)

with Model() as model: # model specifications in PyMC3 are wrapped in a with-statement
    # Define priors
    intercept = Normal('Intercept', 0, sd=20)
    gradient = Normal('gradient', 0, sd=20)


    # Define likelihood
    likelihood = MvNormal('y', mu=intercept + gradient * x,
                        tau=1./(np.stack((y_err, x_err))**2.), observed=y)

    # start the mcmc!
    start = find_MAP() # Find starting value by optimization
    step = NUTS(scaling=start) # Instantiate MCMC sampling algorithm
    trace = sample(2000, step, start=start, progressbar=False) # draw 2000 posterior samples using NUTS sampling

Это поднимает ошибку: LinAlgError: Last 2 dimensions of the array must be square

Так что я пытаюсь пройти MvNormal измеренные значения для х и у (mus) и связанные с ними неопределенности измерений (y_err а также x_err). Но похоже, что это не нравится 2d tau аргумент.

Есть идеи? Это должно быть возможно

Спасибо

Answer 1

Вы можете попробовать адаптировать следующую модель. Является ли "регулярной" линейной регрессией. Но x а также y были заменены гауссовыми распределениями. Здесь я предполагаю не только измеренные значения входных и выходных переменных, но также надежную оценку их погрешности (например, предоставленной измерительным устройством). Если вы не доверяете этим значениям ошибок, вы можете вместо этого попытаться оценить их по данным.

with pm.Model() as model:
    intercept = pm.Normal('intercept', 0, sd=20)
    gradient = pm.Normal('gradient', 0, sd=20)
    epsilon = pm.HalfCauchy('epsilon', 5)
    obs_x = pm.Normal('obs_x', mu=x, sd=x_err, shape=len(x))
    obs_y = pm.Normal('obs_y', mu=y, sd=y_err, shape=len(y))

    likelihood = pm.Normal('y', mu=intercept + gradient * obs_x,
                    sd=epsilon, observed=obs_y)

    trace = pm.sample(2000)

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

df_data = pd.DataFrame(data)
cov = df_data.cov()

with pm.Model() as model:
    intercept = pm.Normal('intercept', 0, sd=20)
    gradient = pm.Normal('gradient', 0, sd=20)
    epsilon = pm.HalfCauchy('epsilon', 5)

    obs_xy = pm.MvNormal('obs_xy', mu=df_data, tau=pm.matrix_inverse(cov), shape=df_data.shape)

    yl = pm.Normal('yl', mu=intercept + gradient * obs_xy[:,0],
                    sd=epsilon, observed=obs_xy[:,1])

mu, sds, elbo = pm.variational.advi(n=20000)
step =  pm.NUTS(scaling=model.dict_to_array(sds), is_cov=True)
trace = pm.sample(1000, step=step, start=mu)

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

Для второй модели я использую ADVI для его инициализации. ADVI может быть хорошим способом инициализации моделей, часто он работает намного лучше, чем find_MAP().

Вы также можете проверить этот репозиторий Дэвидом Хоггом. И в книге " Статистическое переосмысление", где МакЭлриат обсуждает проблему выполнения линейной регрессии, включая ошибки во входных и выходных переменных.