[PYTHON] Essayez d'utiliser un langage de programmation probabiliste (Pyro)

Le langage de programmation probabiliste est utilisé pour la modélisation statistique bayésienne. Je pense que la méthode d'utilisation de Stan, BUGS et JAGS en langage R est courante car il existe de nombreux livres. Les bibliothèques de programmation probabiliste de Python incluent PyMC, PyStan et Edward (maintenant TensorFlow Probability), une extension de TensorFlow. Pyro est sorti en 2018. Pyro est une bibliothèque de programmation probabiliste basée sur PyTorch développée par Uber. Récemment, la société a également publié NumPyro basé sur Jax. Article de Pyro: http://jmlr.org/papers/v20/18-403.html Article de NumPyro: https://arxiv.org/abs/1912.11554

Pyro et NumPyro peuvent facilement implémenter le deep learning bayésien, et je pense qu'ils seront appliqués à l'avenir. (Article mis en œuvre)

Cette fois, je vais essayer les bases du tutoriel officiel de Pyro. J'ai également fait référence à l'article de blog suivant. https://www.hellocybernetics.tech/entry/2019/12/08/033824 https://www.hellocybernetics.tech/entry/2019/12/08/193220

Installer Pyro

On suppose que Python est 3.4 ou supérieur et PyTorch est installé. L'environnement de l'auteur était Ubuntu 16.04, CUDA10, Anaconda, Python3.6, Pytorch 1.2.0.

pip install pyro-ppl

Pyro 1.3.1 est installé. De plus, PyTorch a été automatiquement mis à niveau vers la version 1.4.0.

Mise en œuvre (1): Définition et échantillonnage de la distribution de probabilité

Importez et définissez des semences. la semence est nécessaire pour la reproductibilité de l'échantillonnage. Si vous faites pyro.set_rng_seed, il semble que la valeur de départ dans pytorch soit également décidée.

import torch
import pyro
pyro.set_rng_seed(101)

pyro.distributions peut être similaire à torch.distributions. Cette fois, nous définirons le temps et la température avec la distribution de Bernoulli et la distribution normale.

def weather():
    cloudy = pyro.sample('cloudy', pyro.distributions.Bernoulli(0.3))
    cloudy = 'cloudy' if cloudy.item() == 1.0 else 'sunny'
    mean_temp = {'cloudy': 55.0, 'sunny': 75.0}[cloudy]
    scale_temp = {'cloudy': 10.0, 'sunny': 15.0}[cloudy]
    temp = pyro.sample('temp', pyro.distributions.Normal(mean_temp, scale_temp))
    return cloudy, temp.item()

Empilez-le pour définir les ventes de crème glacée. En d'autres termes, il s'agit d'un modèle de probabilité de la hiérarchie.

def ice_cream_sales():
    cloudy, temp = weather()
    expected_sales = 200. if cloudy == 'sunny' and temp > 80.0 else 50.
    ice_cream = pyro.sample('ice_cream', pyro.distributions.Normal(expected_sales, 10.0))
    return ice_cream

Échantillonnons les ventes de crème glacée et montrons-les dans un histogramme.

import matplotlib.pyplot as plt
N = 1000 #Le nombre d'échantillons
x_list = []
for _ in range(N):
    x = ice_cream_sales()
    x_list.append(x)
plt.hist(x_list)

Une distribution bimodale a été obtenue.

Vectorisation avec la méthode pyro.plate

Si vous utilisez la méthode pyro.plate, il semble que vous puissiez échantillonner sans répéter dans l'instruction for par vectorisation. Cela signifie que l'échantillonnage est effectué en parallèle à partir de la même distribution indépendante. J'ai essayé de réécrire le code jusqu'à présent en utilisant pyro.plate. Définition du modèle

def model(N):
    with pyro.plate("plate", N):
        cloudy = pyro.sample('cloudy', pyro.distributions.Bernoulli(0.3))
        mean_temp = 75.0 - 20.0 * cloudy
        scale_temp = 15.0 - 5.0 * cloudy
        temp = pyro.sample('temp', pyro.distributions.Normal(mean_temp, scale_temp))
        expected_sales = 50.0 + 150.0 * (1 - cloudy) * (temp > 80.0)
        ice_cream = pyro.sample('ice_cream', pyro.distributions.Normal(expected_sales, 10.0))
    return ice_cream

Affichez l'histogramme.

import matplotlib.pyplot as plt
N = 1000 #Le nombre d'échantillons
x_list = model(N)
plt.hist(x_list)

J'ai un histogramme similaire. La légère différence est-elle l'effet de l'échantillonnage en parallèle dû à la vectorisation?

En tant que caractéristique importante de pyro, le nom de la variable est spécifié par le premier argument, tel que pyro.sample ("hoge", dist). Cela semble nécessaire pour définir qu'il s'agit d'une variable stochastique globale. (Je ne suis pas sûr des détails.)

Implémentation (2): inférence bayésienne

Dans la section précédente, nous avons découvert comment générer des données à partir d'un modèle probabiliste. Dans cette section, nous allons essayer de déduire un modèle stochastique à partir des données.

import torch
import pyro
pyro.set_rng_seed(101)
import matplotlib.pyplot as plt

En tant que données, préparez un total de 10 données, 1 valeur x 6 et 0 valeur x 4.

data = []
for _ in range(6):
    data.append(torch.tensor(1.0))
for _ in range(4):
    data.append(torch.tensor(0.0))
data = torch.tensor(data) #Faites-en un tenseur pour la manipulation avec Pyro
plt.hist(data)

L'histogramme ressemblera à celui ci-dessus.

Modélisez une distribution de probabilité pour ces données. Cette fois, considérons la distribution de Bernoulli avec le paramètre "latent_fairness" qui suit la distribution beta.

import pyro.distributions as dist
def model(data):
    alpha0 = torch.tensor(10.0)
    beta0 = torch.tensor(10.0)
    f = pyro.sample("latent_fairness", dist.Beta(alpha0, beta0))
    with pyro.plate("data_plate"):
        pyro.sample("obs", dist.Bernoulli(f), obs=data)

Inférer les paramètres qui correspondent à ce modèle. Il existe plusieurs méthodes d'inférence, mais les variantes Bayes et MCMC sont les principales.

Variante Bayes

Parmi les méthodes d'inférence, Pyro semble avoir les fonctions les plus complètes de la variante Bayes. Ceci est probablement dû au fait que les méthodes automatiques de différenciation et d'optimisation de PyTorch, qui sont à la base de Pyro, sont utilisées pour SVI (Probabilistic Variant Inference) de Variant Bayes.

Tout d'abord, préparez une distribution d'approximation variable. Vous pouvez le définir vous-même, ou vous pouvez le faire automatiquement avec la classe pyro.infer.autoguide.guides.

guide = pyro.infer.autoguide.guides.AutoDelta(model)

Effectue une inférence différentielle.

from pyro.optim import Adam
from pyro.infer import SVI, Trace_ELBO
adam_params = {"lr": 0.0005, "betas": (0.90, 0.999)}
optimizer = Adam(adam_params)
pyro.clear_param_store()
svi = SVI(model, guide, optimizer, loss=Trace_ELBO())
n_steps = 1000 #Nombre d'étapes
losses, f  = [], [] #Pour l'enregistrement pour un tracé ultérieur
for step in range(n_steps):
    loss = svi.step(data)
    losses.append(loss)
    f.append(pyro.param("AutoDelta.latent_fairness").item())

Vous pouvez maintenant en déduire. Tracez la progression.

plt.subplot(2,1,1)
plt.plot(losses)
plt.xlabel("step")
plt.ylabel("loss")
plt.subplot(2,1,2)
plt.plot(f)
plt.xlabel("step")
plt.ylabel("latent_fairness")

Le paramètre "latent_fairness" a convergé vers environ 0,5357. La valeur finale peut être affichée ci-dessous.

for name in pyro.get_param_store():
    print("{}: {}".format(name, pyro.param(name)))

AutoDelta.latent_fairness: 0.535700261592865

Affichez la distribution de probabilité prédite à partir du modèle inféré. Prélevez 1000 échantillons.

from pyro.infer import Predictive
predictive_dist = Predictive(model=model, guide=guide,num_samples=1000, return_sites=["obs"])
plt.hist(predictive_dist.get_samples(None)["obs"].view(1,-1))

La distribution est similaire à l'histogramme des données. (Est-ce légèrement différent parce que le nombre initial de données est aussi petit que 10?)

MCMC MCMC est une technique de force brute et prend du temps.

from pyro.infer import　NUTS, MCMC
nuts_kernel = NUTS(model)
mcmc = MCMC(nuts_kernel, num_samples=5000, warmup_steps=1000)
mcmc.run(data)

Vous pouvez maintenant en déduire. Voyez le résultat.

mcmc.summary()

MCMC peut exprimer la distribution de probabilité prédite en passant un échantillon de paramètres à posterior_samples.

f_mcmc = mcmc.get_samples()["latent_fairness"]
from pyro.infer import Predictive
predictive_dist = Predictive(model=model, posterior_samples={"latent_fairness": f_mcmc}, return_sites=["obs"])
plt.hist(predictive_dist.get_samples(None)["obs"].view(1,-1))

Cela a également une distribution similaire à l'histogramme des données.

Autres méthodes d'inférence

Il existe une estimation MAP et une approximation de Laplace. L'estimation MAP est une forme spéciale d'inférence bayésienne, et comme il s'agit d'une estimation ponctuelle, une distribution de probabilité n'est pas nécessaire. Ce n'est pas une inférence mais une estimation car elle ne demande pas de distribution de probabilité. (Article de référence) *** La variation Bayes effectuée ci-dessus est équivalente à l'estimation MAP car la distribution delta est sélectionnée comme distribution d'approximation de la variation. ** ** L'approximation de Laplace semble possible avec AutoLaplace Approximation de l'autoguide.

Applications

Il semble que vous puissiez faire diverses choses telles que le modèle gaussien mixte, le processus gaussien, l'optimisation bayésienne, le NMF, le filtre de Kalman, le modèle de Markov caché, etc.

À la fin

Je ne fais que commencer avec Pyro, mais d'une manière ou d'une autre, j'ai compris comment l'utiliser. Je pense que c'est une technique merveilleuse de pouvoir mettre en œuvre une théorie qui est difficile avec des formules mathématiques si concises.