Une note d'essayer un simple tutoriel MCMC sur PyMC3

Aperçu

L'un des échantillonneurs MCMC gratuits pouvant être utilisés avec Python est PyMC3. L'autre jour. Une lycéenne assise à côté de moi à Stava a parlé de choses comme "Peut-être que c'est devenu PyMC3 et plus rapide que PyMC2 ..." ou "Stan a des paramètres discrets ...", donc le tutoriel officiel https: En essayant //pymc-devs.github.io/pymc3/getting_started/, j'ai écrit des notes sur certaines erreurs et des expériences supplémentaires. Pour des gens comme "Je veux utiliser des baies hiérarchiques en Python, mais il est difficile d'écrire MCMC avec un scratch complet ..." et "Paramètres de dispersion ...". Un étudiant qui n'est pas très familier avec MCMC l'a fait au milieu de la nuit, donc je pense qu'il y a beaucoup de parties difficiles à lire / petites erreurs / manque de compréhension, mais pardonnez-moi s'il vous plaît.

• Dans MCMC, il est courant de jeter l'échantillon au stade précoce, mais dans l'expérience inférieure, l'échantillon est tracé sans rejet en raison du temps de calcul. Par conséquent, les performances, etc. sont à titre indicatif uniquement.

environnement

Mon environnement PC est OS X10.9.5, Python2.7 +. Vous pouvez saisir ce dont vous avez besoin avec pip. pip install git+https://github.com/pymc-devs/pymc3 De plus, lors de l'installation de PyMC3, veuillez noter que la destination du lien de l'URL peut être différente s'il s'agit d'un petit article ancien. C'est pratique d'avoir des Pandas et Patsy, donc je pense qu'il vaut mieux les mettre avec pip. pip install pandas pip install patsy

Traitement de base

Puisque je vais utiliser une partie du tutoriel pour expliquer où l'erreur s'est produite, je vais résumer les parties minimales nécessaires de la première moitié du tutoriel. Dans cet article, la partie dépendante PyMC3 est explicitement écrite sous la forme de pm. Dans le didacticiel, pm. Est omis, mais sinon, c'est le même code.

Tout d'abord, certains importent.

import numpy as np
%pylab inline
np.random.seed(27)
import pymc3 as pm

Créez un exemple de jeu de données comme suit.

alpha, sigma = 1, 1
beta = [1, 2.5]
size = 100
X1 = np.linspace(0, 1, size)
X2 = np.linspace(0,.2, size)
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma

Définissez le modèle.

#Un modèle qui renvoie Y linéairement sur X1 et X2
#Moyens de définir un modèle nommé model
with pm.Model() as model:

    # Priors for unknown model parameters
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10, shape=2)
    sigma = pm.HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

Exécutez MCMC en utilisant le modèle ci-dessus. (Cela a pris 13,6 secondes avec mon processeur)

from scipy import optimize
with model:
    # obtain starting values via MAP
    start = pm.find_MAP(fmin=optimize.fmin_powell)
    # instantiate sampler
    step = pm.Slice(vars=[sigma])
    # draw 5000 posterior samples
    trace = pm.sample(5000, start=start, step=step)

Quant à l'argument de pm.sample (), le premier 5000 correspond à la taille de l'échantillon que vous souhaitez obtenir. Au démarrage suivant, la valeur initiale est définie. Cette fois, l'estimation MAP est utilisée, mais il est possible de ne pas la spécifier. L'échantillonnage de tranche est spécifié à la dernière étape, mais s'il n'est pas spécifié, il est spécifié. ・ Variable binaire: Métropole binaire ・ Variable discrète: Metropolis ・ Variable continue: NUTS Est attribué par défaut. Dans le modèle simple utilisé ici, il n'y avait pratiquement aucune différence selon la méthode d'échantillonnage. Le Hamiltonian MC et NUTS qui peuvent être utilisés avec Stan peuvent également être utilisés avec PyMC3.

La visualisation est très facile, juste un coup avec traceplot comme suit.

pm.traceplot(trace[4000:])

Les informations récapitulatives peuvent également être visualisées comme suit.

pm.summary(trace[4000:])

Par souci de simplicité, seule la distribution normale a été traitée cette fois, mais bien sûr, diverses autres distributions telles que la distribution gamma, la distribution bêta, la distribution binomiale, la distribution de Poisson, etc. peuvent être traitées.

• Pour la relation avec PyMC2, http://breakbee.hatenablog.jp/entry/2014/08/04/031342 etc. a été utile.

Difficulté

Bien que ce soit un PyMC3 si pratique, il y a quelques problèmes. Essayons le tutoriel.

1. J'obtiens une erreur lorsque j'utilise la méthode ano

Je pense que je viens de voir cette discussion connexe sur Twitter l'autre jour, mais dans ce cas, la cause de l'erreur est probablement la dernière des https://github.com/pymc-devs/pymc3/blob/master/pymc3/model.py ligne

theano.config.compute_test_value = 'raise'

(Référence: https://sites.google.com/site/iwanamidatascience/vol1/support_tokushu). Je ne connais pas les détails du cache, mais après avoir importé PyMC3

import theano
theano.config.compute_test_value = 'ignore'

J'ai pu le résoudre en tapant dans mon cœur. Au fait, il semble que vous puissiez écrire «off» au lieu de «ignore», mais dans le calcul glm etc. scratchpad instance has no attribute 'test_value' Je pense qu'il est prudent de le régler sur "ignorer" pour le moment. À propos, il est également «ignoré» sur la page Factorisation de matrice probabiliste des exemples officiels.

2. Gestion problématique des fonctions et des classes auto-créées

2.1 as_op théorie inutile

Dans la section Déterministique arbitraire du didacticiel, la variable déterministe b est trouvée par la fonction crazy_modulo3 () suivante.

import theano.tensor as T 
from theano.compile.ops import as_op

@as_op(itypes=[T.lscalar], otypes=[T.lscalar])
def crazy_modulo3(value):
    if value > 0: 
        return value % 3
    else :
        return (-value + 1) % 3

with pm.Model() as model_deterministic:
    a = pm.Poisson('a', 1)
    b = crazy_modulo3(a)

La partie @ as_op () est comme un sort pour utiliser une fonction qui calcule une variable de type theano.tensor, et elle est écrite sur la ligne immédiatement avant de définir la fonction pour déclarer le type d'entrée / sortie. Ceci est également écrit dans le didacticiel PyMC3.

Theano needs to know the types of the inputs and outputs of a function, which are specified for as_op by itypes for inputs and otypes for outputs.

** Cependant, cela peut être défini comme suit sans utiliser @ as_op. ** **

import theano.ifelse
def symbolf(x):
    return theano.ifelse.ifelse(T.gt(x,0), x%3, (-x+1)%3)
def crazy_modulo3(value):
    return symbolf(value)

Je suis un peu accro à ça, alors je vais vous l'expliquer.

Par exemple, naïf sans utiliser ifelse

def crazy_modulo3(value):
    if value > 0:
        return value % 3
    else :
        return (-value + 1) % 3

Ensuite, je me fâche si je ne peux pas comparer valeur> 0. Aussi,

def crazy_modulo3(value):
    if T.lt(0,value):
        return value % 3
    else :
        return (-value + 1) % 3

Ensuite, aucune erreur ne se produit, mais lorsque je trace la valeur, T.lt (0, valeur) ne devient pas False (Theano n'a pas de dtype booléen. Au lieu de cela, tous les tenseurs booléens sont représentés dans'int8 ' . --http: //deeplearning.net/software/theano/library/tensor/basic.html). Par conséquent, ifelse est utilisé pour la comparaison forcée. C'est un peu pénible en theano, mais maintenant je peux définir la fonction sans utiliser @ as_op.

2.2 Comment tracer les variables déterministes du tracé

Bien qu'il soit correct de définir la variable déterministe b dans le code ci-dessus, le didacticiel ne décrit pas la méthode d'échantillonnage ultérieure. Essayez-le là-bas

with model_deterministic:
    trace = pm.sample(500)
pm.traceplot(trace)

Si vous essayez, seule la valeur de la variable a sera tracée. J'étais également intéressé par la variable b, donc quand j'ai regardé le contenu de pm., J'ai trouvé quelque chose comme pm.Deterministic (). Le reste est la définition du modèle b = crazy_modulo3(a) À b = pm.Deterministic("b",crazy_modulo3(a)) Si vous le changez en et l'échantillonnez, la valeur de b doit être tracée correctement.

2.3 Évitez d'utiliser as_op pour calculer des variables aléatoires

Dans la section Distributions arbitraires du didacticiel, comment définir la distribution de probabilité que vous souhaitez échantillonner est décrite comme suit.

with pm.Model() as model1:
    alpha = pm.Uniform('intercept', -100, 100)
    # Create custom densities
    beta = pm.DensityDist('beta', lambda value: -1.5 * T.log(1 + value**2), testval=0)
    eps = pm.DensityDist('eps', lambda value: -T.log(T.abs_(value)), testval=1)
    # Create likelihood
    like = pm.Normal('y_est', mu=alpha + beta * X, sd=eps, observed=Y)

Cette méthode définie par l'expression lambda et pm.DensityDist n'est pas un problème. Si vous continuez à échantillonner comme suit, cela fonctionnera correctement.

with model1:
    trace = pm.sample(500)
pm.traceplot(trace)

En tant que code avec une signification similaire, le didacticiel utilise la méthode suivante pour définir une classe (** Je ferai remarquer que ce code sera probablement bientôt incorrect **). Cela n'utilise pas d'expressions lambda, il semble donc avoir l'avantage de faciliter l'écriture de fonctions. À l'exception de la version bêta, ce sont les paramètres que j'ai ajoutés pour qu'ils aient la même forme que le modèle1 ci-dessus.

class Beta(pm.distributions.Continuous):
    def __init__(self, mu, *args, **kwargs):
        super(Beta, self).__init__(*args, **kwargs)
        self.mu = mu
        self.mode = mu

    def logp(self, value):
        mu = self.mu
        return beta_logp(value - mu)

@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
    return -1.5 * np.log(1 + (value)**2)

with pm.Model() as model2:
    beta = Beta('slope', mu=0, testval=0)

    #I add other parameters to follow model1 above
    alpha = pm.Uniform('intercept', -100, 100)
    eps = pm.DensityDist('eps', lambda value: -T.log(T.abs_(value)), testval=1)
    like = pm.Normal('y_est', mu=alpha + beta * X, sd=eps, observed=Y)

L'exemple «@ as_op» apparaît dans le calcul bêta. Quand j'ai utilisé cela pour calculer la variable déterministe plus tôt, l'échantillonnage fonctionnait bien, mais qu'en est-il cette fois-ci?

with model2:
    trace = pm.sample(100)
pm.traceplot(trace)

Si vous faites cela, vous obtiendrez probablement une erreur ↓. AttributeError: 'FromFunctionOp' object has no attribute 'grad'

Apparemment, le message introduit par ** @ as_op n'a pas l'attribut de calcul de gradient grad requis pour l'échantillonnage MCMC **. Était-ce correct parce que le gradient n'est pas calculé avec des variables déterministes? Lorsque j'ai cherché sur Google ici, j'ai trouvé la réponse d'une personne formidable: «Je ne sais pas comment y faire face».

The only way to get the automatic gradient computation is by expressing your density in terms of theano operators. as_op creates a blackbox function for which autodiff will not work so there is no way I know of (except numerical differentiation) to make this work. https://github.com/pymc-devs/pymc3/issues/601

Cependant, voici le résultat de 2.1. En d'autres termes, si @ as_op est mauvais, pourquoi ne pas le définir sans @ as_op en premier lieu? C'est une histoire. Donc,

@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
    return -1.5 * np.log(1 + (value)**2)

Est réécrit comme suit.

def beta_logp(value):
    return -1.5 * T.log(1 + (value)**2)

Comme en 2.1, si vous faites attention aux règles lors du traitement des variables de type theano.tensor (dans ce cas, utilisez T.log au lieu de np.log), vous pouvez calculer en toute sécurité. Les résultats sont presque les mêmes (bien qu'il y ait un peu de flou car il s'agit d'échantillonnage).

** En conclusion, lorsque vous utilisez vos propres fonctions et classes, que vous utilisiez des expressions lambda ou non, vous devez faire attention à la syntaxe de l'eano lors du codage. ** **

3. Vitesse et performances chez glm

Enfin, je parlerai de la vitesse et des performances lors du calcul du modèle linéaire généralisé (glm). Les détails peuvent être trouvés sur https://github.com/pymc-devs/pymc3/issues/544, mais "NUTS est souvent lent par défaut" "Metropolis est rapide" "Selon le type de données et de machine, quand Hamiltonian MC est rapide" Il y a aussi une histoire. " Pour référence, j'ai appliqué les données du didacticiel à glm (régression linéaire ordinaire) dans mon environnement.

import pandas as pd
df = pd.DataFrame({'x1': X1, 'x2': X2, 'y': Y})
with pm.Model() as model:
    pm.glm.glm('y ~ x1+x2', df)
    start = pm.find_MAP(fmin=optimize.fmin_powell)
    step = pm.HamiltonianMC()
    trace = pm.sample(1000, start=start, step=step)
pm.traceplot(trace)

・ HamiltonianMC: Terminé en 1,7 seconde. La convergence est également acceptable. (Figure 1) ・ Metropolis: Terminé en 0,3 seconde. La convergence n'est pas bonne. ・ NOIX: seuls 8 échantillons peuvent être échantillonnés même après 5 minutes, trop tard -Échantillonnage manuscrit sans utiliser la classe pm.glm effectué en 1. de cet article: Terminé en 3,2 secondes. La convergence est inférieure au MC hamiltonien, mais pas si mal. (Figure 2)

Paramètres vrais: ʻalpha = 1, beta [0] = 1, beta [1] = 2.5, sigma = 1` (Figure 1)

(Fig.2: La ligne bleue de beta correspond au coefficient de régression de X1, la ligne verte correspond au coefficient de régression de X2, et sigma correspond à sd)

Le résultat était. À propos, lorsque je coupe NUTS au milieu et que je le trace, il semble que la mise à l'échelle des paramètres soit étrange. Comme solution dans ce cas, sur la page des problèmes ci-dessus

C = pm.approx_hessian(model.test_point)
step = pm.NUTS(scaling=C)

Si c'est le cas, il est écrit que la vitesse sortira. Lorsque j'ai mis à jour la version de Sphinx (si vous ne l'avez pas, installez-la avec pip) et que j'ai tapé les deux lignes ci-dessus, on m'a dit que le mineur de tête n'était pas une valeur positive et ne pouvait pas être calculé. LinAlgError: 2-th leading minor not positive definite

Je suis passé par là parce que cela ne semble pas être une solution, mais il semble y avoir une discussion connexe ici → https://github.com/scikit-learn/scikit-learn/issues/2640

À titre personnel, ** Si vous souhaitez simplement utiliser glm, essayez glm par NUTS en utilisant approx_hessian pour la mise à l'échelle, et si cela ne fonctionne pas, essayez glm par HMC. ** **

En outre, la méthode qui n'utilise pas la classe pm.glm comme décrit dans l'exemple de 1. peut améliorer les performances si vous envisagez comme la définition du puits précédent et si vous souhaitez voir le comportement de l'hyper paramètre ou Quand je veux utiliser un modèle bayésien hiérarchique plus compliqué, j'ai pensé que je devrais adopter cela au lieu d'utiliser pm.glm.

finalement

Je l'ai écrit pendant longtemps avec mon propre mémo, mais j'espère que ce sera une petite référence pour l'endroit qui semble créer une dépendance. Chez @ as_op, il y avait quelque chose d'assez difficile, comme perdre une nuit à essayer de redéfinir le diplômé du côté théano.

Aussi, sur la page officielle, j'ai le sentiment que les exemples sont plus intéressants que le tutoriel (Factorisation de matrice probabiliste, analyse de survie, etc.), donc je vais essayer à nouveau quand j'aurai envie de minuit.

• Pour la définition de grad, reportez-vous à http://deeplearning.net/software/theano/extending/op.html. Quand je pense au T.grad initialement défini comme une fonction qui peut être créée avec la fonction theano., l'entrée semble être une liste, mais si je suis ce tutoriel, je ne peux pas définir une fonction dont l'entrée est une liste et abandonner. fait.

Je l'ajouterai le cas échéant.