pymc3의 다 변수 선형 회귀

Question

나는 최근에 emcee 만 독점적으로 사용한 후 pymc3을 학습하기 시작했으며 일부 개념적 문제가 발생합니다.

제 7 장 Hogg's Fitting a model to data으로 연습하고 있습니다. 이것은 임의의 2 차원 불확실성을 갖는 직선에 대한 mcmc 피트를 포함한다. 나는 이것을 emcee에서 아주 쉽게 수행했지만, pymc이 나에게 몇 가지 문제를주고 있습니다.

본질적으로 다 변수 가우스 확률 (gaussian likelihood)을 사용하는 것으로 귀결됩니다.

여기까지 제가 지금까지 있습니다.

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 x와 y에 대한 측정 값 (mu들) 및 관련 측정 불확도 (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()보다 훨씬 잘 작동합니다.

David Hogg가 repository을 확인할 수도 있습니다. 그리고 책 Statistical Rethinking에서 McElreath는 입력 변수와 출력 변수의 오류를 포함하여 선형 회귀 문제를 논의합니다.