[Python] Estimation bayésienne avec Pyro

Nous publions également le fichier ipynb sur github. https://github.com/isuya0508/exercise_Pyro/blob/master/bayesian_inference_bernulli.ipynb

Qu'est-ce que Pyro?

Pyro est une bibliothèque qui gère les modèles probabilistes avec Pytorch comme back-end. Vous pouvez l'installer à partir de pip.

pip install pyro-ppl 

À ce stade, il est nécessaire d'installer Pytorch à l'avance. Veuillez vous référer à la page officielle pour plus de détails.

Essaye le

Cette fois, considérons l'estimation du paramètre $ p $ à partir des données qui suivent la distribution de Bernoulli $ Ber (p) $. Tout d'abord, importez les modules requis.

from collections import Counter

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

import torch
import torch.distributions.constraints as constraints
import pyro
import pyro.distributions as dist
from pyro.optim import SGD, Adam
from pyro.infer import SVI, Trace_ELBO

%matplotlib inline

pyro.set_rng_seed(0)
pyro.enable_validation(True)

Génération de données

Créez des données avec des nombres aléatoires. À ce stade, le type doit être le tenseur de Pytorch.

obs = stats.bernoulli.rvs(0.7, size=30, random_state=1)
obs = torch.tensor(obs, dtype=torch.float32)
obs

> tensor([1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1.,
          1., 1., 0., 0., 1., 1., 0., 0., 1., 1., 1., 0.])

Vérifiez le nombre de 1 dans les données.

Counter(obs.numpy())

> Counter({1.0: 23, 0.0: 7})

Par conséquent, l'estimation la plus probable pour le paramètre $ p $ est 23/30 $ \ Fallingdotseq 0,77 $. [^ 2] Après cela, nous procéderons à l'estimation bayésienne de $ p $.

Définition de la distribution antérieure

Dans l'estimation bayésienne, la distribution a priori des paramètres est supposée et la distribution a posteriori est obtenue en combinant les données observées.

Pour le paramètre $ p $ de la distribution de Bernoulli $ Ber (p) $, il est courant de supposer une distribution bêta comme distribution antérieure [^ 1] Dans Pyro, la méthode model décrit la distribution antérieure et le modèle des données.

def model(data):
    #En supposant une distribution antérieure
    prior = dist.Beta(1, 1)

    #La modélisation des données
    p = pyro.sample('p', prior)
    for i in range(len(data)):
        pyro.sample(f'obs{i}', dist.Bernoulli(p), obs=data[i])

Cette fois, nous supposons $ Beta (1, 1) $, ce qui est en fait cohérent avec une distribution uniforme.

Définition de la distribution postérieure

Décrivez la distribution postérieure avec la méthode «guide». La distribution postérieure doit être une distribution bêta ainsi que la distribution antérieure. À ce moment, donnez une valeur initiale appropriée comme paramètre de la distribution postérieure.

def guide(data):
    #Définition de la distribution postérieure
    alpha_q = pyro.param('alpha_q', torch.tensor(15), constraint=constraints.positive)
    beta_q = pyro.param('beta_q', torch.tensor(15), constraint=constraints.positive)
    posterior = dist.Beta(alpha_q, beta_q)

    pyro.sample('p', posterior)

Ce paramètre de distribution postérieure $ \ alpha, \ beta $ sera trouvé.

Raccord post-distribué

Quant à la méthode d'estimation des paramètres de la distribution a posteriori, nous utiliserons cette fois l'estimation de variation probabiliste. Il semble que Pyro utilise essentiellement cette méthode. Dans cet exemple de distribution de Bernoulli, il est absurde d'utiliser une méthode d'approximation telle que l'estimation de la variation car la distribution postérieure peut être obtenue analytiquement, mais j'utiliserai cette méthode pour la pratique. (À l'origine, c'est la méthode utilisée lorsque la distribution ne peut pas être calculée analytiquement.)

Pour l'implémentation, générez une instance de Gradient Descent (SGD) et une instance de Probabilistic Variant Estimate (SVI). Lorsque vous créez une instance SVI, vous lui donnez le «modèle» et le «guide» définis ci-dessus.

optimizer = SGD(dict(lr=0.0001, momentum=0.9))
svi = SVI(model, guide, optimizer, loss=Trace_ELBO())

Tout ce que vous avez à faire maintenant est de répéter svi.step (obs) pour mettre à jour les paramètres de distribution postérieure.

NUM_STEPS = 2000
pyro.clear_param_store()

history = {
    'loss': [],
    'alpha': [],
    'beta': []
}

for step in range(1, NUM_STEPS + 1):
    loss = svi.step(obs)
    
    history['loss'].append(loss)
    history['alpha'].append(pyro.param('alpha_q').item())
    history['beta'].append(pyro.param('beta_q').item())
    
    if step % 100 == 0:
        print(f'STEP: {step} LOSS: {loss}')

> 
STEP: 100 LOSS: 17.461310371756554
STEP: 200 LOSS: 18.102468490600586
(Omission)
STEP: 1900 LOSS: 17.97793820500374
STEP: 2000 LOSS: 17.95139753818512

Ici, la perte ajustée et les paramètres de distribution postérieure $ \ alpha et \ beta $ sont enregistrés dans history. Le graphique de ces nombres pour chaque étape ressemble à ceci: (Le code source est omis.)

!

Vérifiez le $ \ alpha, \ beta $ finalement obtenu et vérifiez la valeur attendue et la variance de la distribution postérieure.

infered_alpha = pyro.param('alpha_q').item()
infered_beta = pyro.param('beta_q').item()
posterior = stats.beta(infered_alpha_beta, infered_beta_beta)

print(f'alpha: {infered_alpha}')
print(f'beta: {infered_beta}')
print(f'Expected: {posterior.expect()}')
print(f'Variance: {posterior.var()}')

>
alpha: 25.764650344848633
beta: 7.556574821472168
Expected: 0.7732203787899605
Variance: 0.005109101547631603

Tracons la pré-distribution et la post-distribution. (Le code source est omis.)

Il semble que cela puisse être bien estimé.

Résumé

J'ai essayé d'effectuer une estimation bayésienne simple en utilisant Pyro. Cette fois, c'était une simple estimation bayésienne, mais je pense que la force de Pyro est qu'il peut décrire des modèles complexes de manière flexible et simple. Il existe de nombreux exemples de méthodes probabilistes dans la documentation officielle, et le simple fait de les regarder peut être une expérience d'apprentissage.

[^ 2]: L'estimation la plus probable de $ p $ correspond à la moyenne de l'échantillon. [^ 1]: On sait que la distribution postérieure suit également la distribution beta en faisant de la distribution antérieure de $ p $ une distribution beta. Une telle distribution est appelée distribution a priori conjuguée.