[PYTHON] [Pyro] Modélisation statistique par le langage de programmation probabiliste Pyro ② ~ Modèle de régression simple ~

introduction

Dans l'article précédent (https://qiita.com/WKudo/items/efd0841959822572fac1), j'ai présenté un aperçu du langage de programmation probabiliste Pyro et comment échantillonner des échantillons à l'aide de Pyro. Cette fois, en utilisant Pyro, "Introduction à l'analyse des données par modélisation statistique bayésienne à partir de Practical Data Science Series R et Stan" (ci-après dénommé livre de référence) Sur la base du premier exemple (chapitre 3, partie 2, modèle de régression simple) dans (Masu.), Nous allons l'implémenter dans Pyro. L'implémentation est disponible sur 3-2 Simple Regression Model (Google Colaboratory).

Données d'utilisation et exemple d'introduction

Les données traitées cette fois (voir livre 3-2) sont ** Relation entre les ventes de bière et la température **. La modélisation statistique construit un modèle qui explique les ventes en fonction de la température. Des exemples de données réelles et de diagrammes de dispersion sont les suivants.

#Lecture des données
beer_sales_2 = pd.read_csv(
    'https://raw.githubusercontent.com/logics-of-blue/book-r-stan-bayesian-model-intro/master/book-data/3-2-1-beer-sales-2.csv'
) # shape = (100, 2)
beer_sales_2.head()

#Visualisation
beer_sales_2.plot.scatter('temperature', 'sales')
plt.title("Relation entre les ventes de bière et la température")

Il semble y avoir une corrélation positive en quelque sorte. Ce qui suit est une description d'un modèle qui explique la vente en termes de température à l'aide d'un modèle de régression simple. sales_i \sim \rm{Normal}(Intercept + \beta * temperature_i, \sigma^2) Dans la modélisation statistique bayésienne, les paramètres $ Intercept $, $ \ beta $, $ \ sigma ^ 2 $ que vous souhaitez estimer dans l'équation ci-dessus sont traités comme des variables stochastiques (appelées variables latentes), et les données d'observation sont utilisées après coup. Estimez la distribution.

Implémentation par Pyro

Dans Pyro, la description du modèle et l'estimation des paramètres sont effectuées selon la procédure suivante.

1. description de la méthode modèle
2. Estimation de la distribution postérieure Dans ce qui suit, nous expliquerons les composants un par un.

description de la méthode du modèle

Dans Pyro, le processus par lequel le modèle statistique supposé génère des données est décrit dans une méthode. Le nom de la méthode est généralement «model». À ce moment-là, veillez à respecter les points suivants.

• Décrivez le processus de génération de données à l'aide de «pyro.sample» et de «pyro.plate» présenté dans Article précédent.
• Reçoit les données d'observation comme argument et décrit où correspondent les données d'observation dans le processus de génération de données. Voici la mise en œuvre.

import torch
import pyro
import pyro.distributions as dist


def model(temperature: torch.Tensor, sales: torch.Tensor=None):
    intercept = pyro.sample("intercept", dist.Normal(0, 1000))
    beta = pyro.sample("beta", dist.Normal(0, 1000))
    sigma2 = pyro.sample("sigma2", dist.Uniform(0, 1000))
    with pyro.plate("plate", size=temperature.size(0)):
        y_ = intercept + beta * temperature
        y = pyro.sample("y", dist.Normal(y_, sigma2), obs=sales)
    return y

Les quatre variables ($ Intercept $, $ \ beta $, $ \ sigma ^ 2 $, $ y $) définies par pyro.sample dans model sont traitées comme des variables stochastiques. Ensuite, parmi les variables de probabilité, celles qui ne correspondent pas aux données observées (autres que $ y $) sont des variables latentes, et seront la cible pour estimer ultérieurement la distribution a posteriori.

Estimation de la distribution postérieure

Après avoir écrit la méthode du modèle pour le modèle statistique, la distribution a posteriori de la variable latente est estimée. Les deux méthodes d'estimation suivantes sont principalement fournies.

1. MCMC
2. Raisonnement variable

MCMC Les méthodes requises pour l'estimation sont fournies dans «pyro.infer». Pour faire une estimation à l'aide du noyau NUTS, écrivez comme suit.

import pyro.infer as infer

#Variables explicatives / variables observées torche.Convertir en tenseur
temperature = torch.Tensor(beer_sales_2.temperature)
sales = torch.Tensor(beer_sales_2.sales)

#Estimation de la distribution postérieure
nuts_kernel = infer.NUTS(model, adapt_step_size=True, jit_compile=True, ignore_jit_warnings=True)
mcmc = infer.MCMC(nuts_kernel, num_samples=3000, warmup_steps=200, num_chains=3)
mcmc.run(temperature, sales)

Parmi les variables stochastiques de la méthode du modèle, la distribution a posteriori de la section (intersection), de la pente (bêta) et de la variance (sigma2) à estimer a été estimée. La section est d'environ 21 et la pente est d'environ 2,5, ce qui est presque le même que le résultat indiqué dans le livre de référence.

Estimation par inférence variable

Pyro peut également estimer par inférence variable. Dans l'inférence de variable, la variable latente que vous voulez trouverZDistribution postérieure dep(Z|X)Vers une autre fonctionq(Z)Pensez à vous rapprocher de. Dans Pyro ceciq(Z)Méthode correspondant à(guide)Si vous écrivezp(Z|X)Meilleur approximatifq(Z)Peut être calculé. Organisons la situation en la remplaçant par un exemple. En regardant en arrière le processus de génération de données décrit par la méthode model, nous pouvons voir qu'il y a trois variables latentes que nous voulons estimer pour ce problème: ʻintercept, betaetsigma2. Le processus de génération de ces trois variables latentes peut être décrit dans la méthode «guide». Cependant, si vous rédigez le guide de manière appropriée$q(Z)$À$p(Z|X)$Je ne peux pas m'en approcher. Dans Pyropyro.param`Les paramètres déclarés parq(Z)Quandp(Z|X)L'inférence de variante est réalisée en mettant à jour la distribution de. Pour ce faire, nous utilisons la différenciation automatique et la descente de gradient de Pytorch. Sur la base de ce qui précède, comment écrire la méthode «guide» est résumée ci-dessous.

• Déclarez le paramètre cible de mise à jour et la valeur initiale avec pyro.param.
• Lisez le paramètre dans la méthode guide. Les contraintes sont spécifiées par dist.constraints.
• Générer des variables latentes avec dist. Distribution de probabilité (paramètre cible de mise à jour) et pyro.sample.
• L'argument doit être le même que la méthode model.

Voici un exemple de l'implémentation de «guide» pour ce problème. Cette fois, nous avons supposé que les trois variables latentes étaient indépendantes et avons choisi la distribution delta comme forme fonctionnelle.

#Déclaration des variables et réglage des valeurs initiales
pyro.param("intercept_q", torch.tensor(0.))
pyro.param("beta_q", torch.tensor(0.))
pyro.param("sigma2_q", torch.tensor(10.))

# q(Z)Implémentation de
def guide(temperature, sales=None):
    intercept_q = pyro.param("intercept_q")
    beta_q = pyro.param("beta_q")
    sigma2_q = pyro.param("sigma2_q", constraint=dist.constraints.positive) 
    intercept = pyro.sample("intercept", dist.Delta(intercept_q))
    beta = pyro.sample("beta", dist.Delta(beta_q))
    sigma2 = pyro.sample("sigma2", dist.Delta(sigma2_q))

Par exemple, si vous voulez une distribution normale au lieu d'une distribution delta, définissez la moyenne ($ \ mu ) et la variance ( \ sigma ^ 2 ) dans `pyro.param` pour chaque variable latente, puis` dist.Normal. Tout ce que vous avez à faire est de le réécrire pour qu'il soit généré à partir de `( \ mu $, $ \ sigma ^ 2 $) en utilisant pyro.sample. De cette manière, vous pouvez modifier de manière flexible la méthode d'estimation de la distribution postérieure en modifiant la pièce de guidage. À propos, si toutes les variables latentes ont une distribution delta indépendante comme dans l'exemple de code, guide = infer.autoguide.guides.AutoDelta(model) C'est la même chose même si vous écrivez. ʻInfer.autogide.guides fournit une fonction pratique qui crée automatiquement une méthode guide`. (http://docs.pyro.ai/en/0.2.1-release/contrib.autoguide.html)

Si vous écrivez les méthodes model et guide, vous pouvez facilement trouver la distribution postérieure.

#Une méthode pour conserver la valeur du paramètre à estimer pour chaque époque
def update_param_dict(param_dict):
    for name in pyro.get_param_store():
        param_dict[name].append(pyro.param(name).item())

#Exécution de l'inférence de variantes
adam_params = {"lr": 1e-1, "betas": (0.95, 0.999)}
oprimizer = pyro.optim.Adam(adam_params)

svi = infer.SVI(model, guide, oprimizer, loss=infer.JitTrace_ELBO(),)

n_steps = 10000
losses = []
param_dict = defaultdict(list)

for step in tqdm(range(n_steps)):
    loss = svi.step(temperature, sales,)
    update_param_dict(param_dict)
    losses.append(loss)

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

    intercept_q: 21.170936584472656
    beta_q: 2.459129810333252
    sigma2_q: 16.684457778930664

Avec environ 21 sections, environ 2,5 pentes et environ 17 dispersions, des résultats d'estimation raisonnables ont été obtenus. D'après la figure ci-dessous, on peut voir que la perte et chaque paramètre cible de mise à jour sont également convergés.

Accélération par GPU

Maintenant que nous avons trouvé que la distribution postérieure peut être estimée par inférence de variables, nous aimerions considérer l'évolutivité en utilisant le GPU, qui est le plus grand avantage de Pyro. Le moyen le plus simple d'utiliser le GPU avec Pyro est au début du code torch.set_default_tensor_type(torch.cuda.FloatTensor) En écrivant, tous les Tensors générés sont placés sur le GPU depuis le début. À propos, dans l'exemple, le nombre d'échantillons était de 100, mais si cela augmente à 1000, 10000, ..., comment le temps nécessaire à l'inférence de variables augmentera-t-il? Ici, en supposant que les sections, les pentes et les variances estimées à partir de 100 échantillons sont des valeurs vraies, 1000, 10000, ... échantillons sont générés à partir de la distribution vraie (provisoire), et les échantillons sont générés. Nous avons comparé le temps requis pour l'inférence de variables avec le CPU et le GPU comme données d'observation.

Comme vous pouvez le voir dans le tableau, dans la zone où le nombre d'échantillons est de 10 000 ou moins, l'avantage de la parallélisation n'est pas reçu et le processeur est plus rapide. Cependant, la différence a commencé à s'ouvrir lorsque le nombre d'échantillons dépassait 100 000, et le résultat était qu'il y avait une différence de 30 fois ou plus pour des pièces de 10 $ ^ 8 $. On peut voir que la modélisation statistique peut être effectuée sans stress, même lorsqu'il s'agit de données extrêmement volumineuses à l'aide de GPU.

Résumé

Cette fois, j'ai introduit une méthode d'analyse par un modèle de régression simple utilisant Pyro. Le processus par lequel le modèle statistique supposé génère des données peut être décrit dans la méthode du «modèle», et la distribution postérieure peut être estimée par MCMC ou par inférence différentielle. Nous avons également mené des expériences sur l'accélération de l'utilisation de GPU et confirmé que Pyro peut effectuer une analyse évolutive.