[PYTHON] [Pyro] Statistische Modellierung mit der probabilistischen Programmiersprache Pyro ~ Verteiltes Analysemodell, normales lineares Modell ~

Einführung

Zurück ・ Letztes Mal, die Grundlagen der probabilistischen Programmiersprache Pyro Ich stellte vor, wie man es benutzt und wie man mit Pyro die Parameter eines einfachen Regressionsmodells ableitet. Diesmal das "Distributed Analysis Model", das in Kapitel 6 von Teil 3 von Reference Books und das "Normal Linear Model", das in Kapitel 7 erscheint. Ich werde vorstellen, wie man mit Pyro "" realisiert.

Verteiltes Analysemodell

Nutzungsdaten & Beispieleinführung

Ähnlich wie in Vorheriger Artikel beziehen sich die diesmal verarbeiteten Daten auch auf den Bierverkauf. Das letzte Mal war die erklärende Variable die Temperatur, aber dieses Mal ist die erklärende Variable ** Wetter **.

Das ** verteilte Analysemodell ** wird als Modell zur Erklärung des Bierverkaufs verwendet, wobei die qualitative Wettervariable als erklärende Variable verwendet wird. Im Nachschlagewerk wird es wie folgt erklärt.

Bei diesem Modell wird davon ausgegangen, dass sich der Durchschnittswert der Antwortvariablen je nach Wetterlage ändert. Ein Modell, das qualitative Daten als erklärende Variablen verwendet und eine Normalverteilung für die Wahrscheinlichkeitsverteilung annimmt, wird als Dispersionsanalysemodell bezeichnet.

(Shinya Baba. Praktische Data Science-Reihe, beginnend mit R und Stand Einführung in die Datenanalyse durch Bayesianische statistische Modellierung (japanische Ausgabe) (Kindle-Position Nr. 3608-3610). Kindle-Ausgabe.)

In diesem Problem gibt es drei Arten von Wetter: {sonnig, regnerisch, bewölkt}, sodass eine Dummy-Variable ($ \ it {sonnig} ) für sonnig und eine Dummy-Variable ( \ it {regnerisch} ) für Regen enthalten sind. ， Das Modell ist wie folgt. Schätzen Sie die latenten Variablen ( 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)

Implementierung

Das Folgende zeigt eine Implementierung, die die posteriore Verteilung durch MCMC schätzt.

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

#Daten gelesen
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')
#Dummy-Variable
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)

#Schätzung der posterioren Verteilung durch 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()

Sie können sehen, dass es mit den Ergebnissen im Buch übereinstimmt. Da das 95% -Konfidenzintervall von $ beta_1 $ zwischen 14,12 und 25,27 liegt, können wir sehen, dass der Umsatz bei sonnigem Wetter tendenziell steigt, wenn es bewölkt ist. Andererseits kann, da das Konfidenzintervall von $ beta_2 $ 0 überspannt, gesehen werden, dass es keinen großen Unterschied gibt, wenn es regnet als wenn es bewölkt ist.

Normales lineares Modell

Nutzungsdaten & Beispieleinführung

Das einfache Regressionsmodell geht von einem Modell aus, das den Umsatz nach Temperatur erklärt, und das Dispersionsanalysemodell geht von einem Modell aus, das den Umsatz nach Wetterlage erklärt. Aus den bisherigen Modellierungsergebnissen ging hervor, dass sowohl Temperatur als auch Wetter umsatzbezogene Faktoren sind. Um den Umsatz zu erklären, betrachten wir hier die Erstellung eines Modells, das sowohl die quantitative Temperaturvariable als auch die qualitative Wettervariable berücksichtigt. In dieser Situation wird das ** normale lineare Modell ** verwendet. Die Erklärung durch das Nachschlagewerk lautet wie folgt.

Unabhängig von quantitativen oder qualitativen Daten wird ein Modell, in dem mehrere erklärende Variablen für einen linearen Prädiktor verwendet werden, eine Gleichheitsfunktion eine Verknüpfungsfunktion und eine Normalverteilung für eine Wahrscheinlichkeitsverteilung im Allgemeinen als normales lineares Modell bezeichnet. ..

Das Modell wird wie folgt beschrieben. Schätzen Sie die latenten Variablen ($ 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)

Implementierung

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

#Daten gelesen
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')
#Datentransformation
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)
#Schätzung durch variable Inferenz
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)

#Anzeige der Schätzergebnisse
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

Der Abschnitt ist ungefähr 20,2, $ beta_1 $ ist ungefähr 29,5, $ beta_2 $ ist ungefähr -3,54 und $ beta_3 $ ist ungefähr 2,55, was mit den MCMC-Schätzergebnissen im Nachschlagewerk übereinstimmt.

Zusammenfassung

Zurück und in diesem Artikel ist Pyro ein lineares Regressionsmodell (einzelnes Regressionsmodell, Dispersionsanalysemodell, normales lineares Modell), das eine Normalverteilung für Antwortvariablen annimmt. Ich habe vorgestellt, wie man es realisiert. Ab dem nächsten Mal werde ich eine Methode zur Realisierung der Modellierung durch ein Poisson-Regressionsmodell unter der Annahme einer Poisson-Verteilung für die Verteilung von Antwortvariablen mit Pyro vorstellen.