[PYTHON] [Pyro] Modélisation statistique par le langage de programmation probabiliste Pyro ③ ~ Modèle d'analyse distribuée, modèle linéaire normal ~

introduction

Précédent ・ Précédent, les bases du langage de programmation probabiliste Pyro J'ai présenté comment l'utiliser et comment inférer bayésien les paramètres d'un modèle de régression simple avec Pyro. Cette fois, le «modèle d'analyse distribuée» qui apparaît au chapitre 6 de la partie 3 de Ouvrages de référence et le «modèle linéaire normal» qui apparaît au chapitre 7. Je vais vous présenter comment réaliser "" avec Pyro.

Modèle d'analyse distribué

Données d'utilisation et exemple d'introduction

Similaire à Article précédent, les données traitées cette fois-ci sont également liées aux ventes de bière. La dernière fois, la variable explicative était la température, mais cette fois la variable explicative est ** météo **.

Le ** modèle d'analyse distribuée ** est utilisé comme modèle pour expliquer les ventes de bière, en utilisant la variable qualitative du temps comme variable explicative. Dans le livre de référence, il est expliqué comme suit.

Ce modèle suppose que la valeur moyenne des variables de réponse change en fonction de la météo. Un modèle qui utilise des données qualitatives comme variables explicatives et suppose une distribution normale pour la distribution de probabilité est appelé un modèle d'analyse de variance.

(Shinya Baba. Série pratique sur la science des données commençant par R et Stand Introduction à l'analyse des données par modélisation statistique bayésienne (édition japonaise) (Kindle Position n ° 3608-3610).

Dans ce problème, il existe trois types de temps, {ensoleillé, pluvieux, nuageux}, donc une variable fictive ($ \ it {ensoleillé} ) représentant ensoleillé et une variable fictive ( \ it {pluvieux} ) représentant la pluie sont incorporées. ， Le modèle est le suivant. Estimez les variables latentes ( Intercept, \ beta_1, \ beta_2, sigma ^ 2 $).

sales_i \sim \rm{Normal}(Intercept + \beta_1\it{sunny_i} + \beta_2\it{rainy}_i , sigma^2)

la mise en oeuvre

Ce qui suit montre une implémentation qui estime la distribution postérieure par MCMC.

import pyro
import pyro.distributions as dist
import pyro.infer as infer


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

#Lecture des données
beer_sales_3 = pd.read_csv('https://github.com/logics-of-blue/book-r-stan-bayesian-model-intro/raw/master/book-data/3-6-1-beer-sales-3.csv')
#Variable muette
rainy_sunny = pd.get_dummies(beer_sales_3.weather, drop_first=True)
sunny = torch.Tensor(rainy_sunny.sunny)
rainy = torch.Tensor(rainy_sunny.rainy)
sales = torch.Tensor(beer_sales_3.sales)

#Estimation de la distribution postérieure par MCMC
nuts_kernel = infer.NUTS(model, adapt_step_size=True, jit_compile=True, ignore_jit_warnings=True)
mcmc = infer.MCMC(nuts_kernel, num_samples=5000, warmup_steps=200, num_chains=1)
mcmc.run(sunny, rainy, sales)
mcmc.summary()

Vous pouvez voir que cela est cohérent avec les résultats du livre. Puisque l'intervalle de confiance à 95% de $ beta_1 $ est de 14,12 à 25,27, nous pouvons voir que les ventes ont tendance à augmenter lorsque le temps est ensoleillé par rapport à lorsqu'il est nuageux. D'un autre côté, puisque l'intervalle de confiance de $ beta_2 $ chevauche 0, on peut voir qu'il n'y a pas de grande différence quand il pleut que quand il fait nuageux.

Modèle linéaire normal

Données d'utilisation et exemple d'introduction

Dans le modèle de régression simple, nous avons supposé un modèle qui explique les ventes en fonction de la température, et dans le modèle d'analyse de dispersion, nous avons supposé un modèle qui explique les ventes en fonction des conditions météorologiques. D'après les résultats de la modélisation jusqu'à présent, il a été constaté que la température et le temps sont des facteurs liés aux ventes. Ici, pour expliquer les ventes, nous envisageons de construire un modèle qui intègre à la fois la variable quantitative de la température et la variable qualitative du temps. Le ** modèle linéaire normal ** est utilisé dans cette situation. L'explication par le livre de référence est la suivante.

Indépendamment des données quantitatives ou qualitatives, un modèle dans lequel plusieurs variables explicatives sont utilisées pour un prédicteur linéaire, une fonction égale est une fonction de lien et une distribution normale est utilisée pour une distribution de probabilité est généralement appelé un modèle linéaire normal. ..

Le modèle est décrit comme suit. Estimez les variables latentes ($ Intercept, \ beta_1, \ beta_2, \ beta_3, sigma ^ 2 $).

sales_i \sim \rm{Normal}(Intercept + \beta_1\it{sunny_i} + \beta_2\it{rainy_i} + \beta_3\it{temperature_i}, sigma^2)

la mise en oeuvre

import pyro
import pyro.distributions as dist
import pyro.infer as infer


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

#Lecture des données
beer_sales_4 = pd.read_csv('https://github.com/logics-of-blue/book-r-stan-bayesian-model-intro/raw/master/book-data/3-7-1-beer-sales-4.csv')
#Transformation des données
rainy_sunny = pd.get_dummies(beer_sales_4.weather, drop_first=True)
sunny = torch.Tensor(rainy_sunny.sunny)
rainy = torch.Tensor(rainy_sunny.rainy)
temperature = torch.Tensor(beer_sales_4.temperature)
sales = torch.Tensor(beer_sales_4.sales)
#Estimation par inférence variable
guide = infer.autoguide.guides.AutoDelta(model)
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 = 100000
losses = []
param_dict = defaultdict(list)
for step in tqdm(range(n_steps)):
    loss = svi.step(sunny, rainy, temperature ,sales,)
    update_param_dict(param_dict)
    losses.append(loss)

#Affichage des résultats d'estimation
for name in pyro.get_param_store():
    print("{}: {}".format(name, pyro.param(name)))

AutoDelta.intercept: 20.227853775024414 AutoDelta.beta1: 29.45983123779297 AutoDelta.beta2: -3.535917043685913 AutoDelta.beta3: 2.547081708908081 AutoDelta.sigma2: 15.715911865234375

La section est d'environ 20,2, $ beta_1 $ est d'environ 29,5, $ beta_2 $ est d'environ -3,54 et $ beta_3 $ est d'environ 2,55, ce qui est cohérent avec les résultats d'estimation MCMC dans le livre de référence.

Résumé

Précédent et dans cet article, Pyro est un modèle de régression linéaire (modèle de régression unique, modèle d'analyse de dispersion, modèle linéaire normal) qui suppose une distribution normale pour les variables de réponse. J'ai présenté comment le réaliser. À partir de la prochaine fois, j'introduirai une méthode pour réaliser la modélisation par modèle de régression de Poisson en supposant une distribution de Poisson pour la distribution des variables de réponse avec Pyro.