Régression linéaire multivariée dans pymc3

Question

J'ai récemment commencé à apprendre pymc3 après avoir utilisé exclusivement emcee pour les âges et je rencontre des problèmes conceptuels. Je m'entraîne avec le chapitre 7 de Hogg's Fitting a model to data. Ceci implique un ajustement de mcmc à une ligne droite avec des incertitudes 2d arbitraires. Je l'ai fait assez facilement dans emcee, mais pymc me donne quelques problèmes.

Il se résume essentiellement à l'utilisation d'une probabilité gaussienne multivariée.

Voici ce que j'ai jusqu'à présent.

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 

Cela soulève l'erreur: LinAlgError: Last 2 dimensions of the array must be square

Je suis en train de passer MvNormal les valeurs mesurées pour x et y (mu s) et leurs incertitudes de mesure associées (y_err et x_err). Mais il semble qu'il n'aime pas l'argument 2d tau.

Des idées? Cela doit être possible

Merci

Answer 1

Vous pouvez essayer en adaptant le modèle suivant. Est une régression linéaire "régulière". Mais x et y ont été remplacés par des distributions gaussiennes. Ici, je suppose non seulement les valeurs mesurées des variables d'entrée et de sortie, mais aussi une estimation fiable de leur erreur (par exemple fournie par un appareil de mesure). Si vous ne faites pas confiance à ces valeurs d'erreur, vous pouvez essayer de les estimer à partir des données.

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) 

Si vous estimer l'erreur à partir des données qu'il pourrait être raisonnable de supposer qu'ils pourraient être mis en corrélation et, par conséquent, au lieu d'utiliser deux gaussienne, vous pouvez utiliser une gaussienne multivariée séparée. Dans ce cas, vous finirez avec un modèle comme ce qui suit:

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) 

Notez que dans le modèle précédent, la matrice de covariance a été calculée à partir des données. Si vous allez faire cela, je pense qu'il vaut mieux utiliser le premier modèle, mais si vous allez plutôt estimer la matrice de covariance, alors le deuxième modèle pourrait être une approche sensée.

Pour le second modèle, j'utilise ADVI pour l'initialiser. ADVI peut être un bon moyen d'initialiser les modèles, souvent cela fonctionne beaucoup mieux que find_MAP().

Vous pouvez également vérifier ceci repository par David Hogg. Et le livre Statistical Rethinking où McElreath discuter du problème de faire la régression linéaire, y compris les erreurs dans les variables d'entrée et de sortie.