[PYTHON] [Pyro] Statistische Modellierung mit der probabilistischen Programmiersprache Pyro ~ Einfaches Regressionsmodell ~

Einführung

Im vorherigen Artikel (https://qiita.com/WKudo/items/efd0841959822572fac1) habe ich einen Überblick über die probabilistische Programmiersprache Pyro und das Abtasten von Samples mit Pyro gegeben. Diesmal unter Verwendung von Pyro "Einführung in die Datenanalyse durch statistische Bayes'sche Modellierung, beginnend mit den praktischen Datenwissenschaftsreihen R und Stan" (im Folgenden als Nachschlagewerk bezeichnet) Basierend auf dem ersten Beispiel (Kapitel 3, Teil 2, Einfaches Regressionsmodell) in (Masu.) Werden wir es in Pyro implementieren. Die Implementierung ist unter 3-2 Simple Regression Model (Google Colaboratory) verfügbar.

Nutzungsdaten & Beispieleinführung

Die diesmal verarbeiteten Daten (siehe Buch 3-2) sind ** Beziehung zwischen Bierverkauf und Temperatur **. Die statistische Modellierung erstellt ein Modell, das den Umsatz anhand der Temperatur erklärt. Beispiele für tatsächliche Daten und Streudiagramme sind wie folgt.

#Daten gelesen
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()

#Visualisierung
beer_sales_2.plot.scatter('temperature', 'sales')
plt.title("Beziehung zwischen Bierverkauf und Temperatur")

Es scheint irgendwie eine positive Korrelation zu geben. Das Folgende ist eine Beschreibung eines Modells, das den Verkauf anhand eines einfachen Regressionsmodells in Bezug auf die Temperatur erklärt. sales_i \sim \rm{Normal}(Intercept + \beta * temperature_i, \sigma^2) Bei der Bayes'schen statistischen Modellierung werden die in der obigen Gleichung zu schätzenden Parameter $ Intercept $, $ \ beta $, $ \ sigma ^ 2 $ als stochastische Variablen behandelt (diese werden als latente Variablen bezeichnet), und die Beobachtungsdaten werden nachträglich verwendet. Schätzen Sie die Verteilung.

Implementierung durch Pyro

In Pyro werden die Modellbeschreibung und die Parameterschätzung gemäß dem folgenden Verfahren durchgeführt.

1. Beschreibung der Modellmethode
2. Schätzung der posterioren Verteilung Im Folgenden werden die Komponenten einzeln erläutert.

Beschreibung der Modellmethode

In Pyro wird der Prozess, mit dem das angenommene statistische Modell Daten generiert, in einer Methode beschrieben. Der Name der Methode lautet normalerweise "model". Achten Sie zu diesem Zeitpunkt darauf, die folgenden Punkte zu erfüllen.

• Beschreiben Sie den Datengenerierungsprozess mit "pyro.sample" und "pyro.plate", die in Vorheriger Artikel eingeführt wurden.
• Empfängt Beobachtungsdaten als Argument und beschreibt, wo die Beobachtungsdaten im Datengenerierungsprozess übereinstimmen. Unten ist die Implementierung.

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

Alle vier Variablen ($ Intercept $, $ \ beta $, $ \ sigma ^ 2 $, $ y $), die durch pyro.sample in model definiert sind, werden als stochastische Variablen behandelt. Dann sind unter den Wahrscheinlichkeitsvariablen diejenigen, die nicht den beobachteten Daten entsprechen (außer $ y $), latente Variablen und werden das Ziel für die spätere Schätzung der posterioren Verteilung sein.

Schätzung der posterioren Verteilung

Nach dem Schreiben der Modellmethode für das statistische Modell wird die posteriore Verteilung der latenten Variablen geschätzt. Die folgenden zwei Schätzmethoden werden hauptsächlich bereitgestellt.

1. MCMC
2. Argumentationsvariante

MCMC Die für die Schätzung erforderlichen Methoden sind in "pyro.infer" angegeben. Um eine Schätzung mit dem NUTS-Kernel vorzunehmen, schreiben Sie wie folgt.

import pyro.infer as infer

#Erklärende Variablen / beobachtete Variablen Fackel.In Tensor konvertieren
temperature = torch.Tensor(beer_sales_2.temperature)
sales = torch.Tensor(beer_sales_2.sales)

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

Unter den stochastischen Variablen in der Modellmethode wurde die posteriore Verteilung des zu schätzenden Abschnitts (Achsenabschnitt), der Steigung (Beta) und der Varianz (Sigma2) geschätzt. Der Abschnitt ist ungefähr 21 und die Steigung ist ungefähr 2,5, was fast dem Ergebnis entspricht, das im Nachschlagewerk gezeigt wird.

Schätzung durch variable Inferenz

Pyro kann auch durch variable Inferenz schätzen. In der Variableninferenz die latente Variable, die Sie suchen möchtenZPosteriore Verteilung vonp(Z|X)Zu einer anderen Funktionq(Z)Betrachten Sie die Annäherung an. In Pyro dasq(Z)Methode entsprechend(guide)Wenn du schreibstp(Z|X)Beste ungefähreq(Z)Kann berechnet werden. Lassen Sie uns die Situation organisieren und durch ein Beispiel ersetzen. Wenn wir auf den Datengenerierungsprozess zurückblicken, der durch die "Modell" -Methode beschrieben wird, können wir sehen, dass es drei latente Variablen gibt, die wir für dieses Problem schätzen möchten: "Intercept", "Beta" und "Sigma2". Der Prozess der Erzeugung dieser drei latenten Variablen kann in der "Guide" -Methode beschrieben werden. Wenn Sie den Leitfaden jedoch entsprechend schreibenq(Z)Zup(Z|X)Ich kann nicht näher kommen. In Pyropyro.paramDie von deklarierten Parameterq(Z)Wannp(Z|X)Varianteninferenz wird durch Aktualisieren der Verteilung von realisiert. Dabei verwenden wir die automatische Differenzierung und den Gradientenabstieg von Pytorch. Auf der Grundlage der obigen Ausführungen wird im Folgenden zusammengefasst, wie die "Guide" -Methode geschrieben wird.

• Deklarieren Sie den Aktualisierungszielparameter und den Anfangswert mit pyro.param.
• Lesen Sie den Parameter in der guide Methode. Einschränkungen werden durch "dist.constraints" angegeben.
• Generieren Sie latente Variablen mit "dist. Wahrscheinlichkeitsverteilung (Zielparameter aktualisieren)" und "pyro.sample".
• Das Argument muss mit der Methode "model" identisch sein.

Das Folgende ist ein Beispiel für die Implementierung von "Anleitung" für dieses Problem. Dieses Mal nahmen wir an, dass die drei latenten Variablen unabhängig waren, und wählten die Deltaverteilung als funktionale Form.

#Deklaration von Variablen und Einstellung von Anfangswerten
pyro.param("intercept_q", torch.tensor(0.))
pyro.param("beta_q", torch.tensor(0.))
pyro.param("sigma2_q", torch.tensor(10.))

# q(Z)Implementierung von
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))

Wenn Sie beispielsweise eine Normalverteilung anstelle einer Deltaverteilung wünschen, definieren Sie den Mittelwert ($ \ mu ) und die Varianz ( \ sigma ^ 2 ) in `pyro.param` für jede latente Variable und dann` dist.Normal. Alles was Sie tun müssen, ist es neu zu schreiben, so dass es aus `( \ mu $, $ \ sigma ^ 2 $) mit pyro.sample generiert wird. Auf diese Weise können Sie die Methode zur Schätzung der posterioren Verteilung flexibel ändern, indem Sie den Führungsteil ändern. Übrigens, wenn alle latenten Variablen eine unabhängige Deltaverteilung haben, wie im Codebeispiel, guide = infer.autoguide.guides.AutoDelta(model) Es ist das gleiche, auch wenn Sie schreiben. infer.autogide.guides bietet eine praktische Funktion, die automatisch eine guide -Methode erstellt. (http://docs.pyro.ai/en/0.2.1-release/contrib.autoguide.html)

Wenn Sie die Methoden "Modell" und "Leitfaden" schreiben, können Sie die hintere Verteilung leicht finden.

#Eine Methode zum Halten des Werts des Parameters, der für jede Epoche geschätzt werden soll
def update_param_dict(param_dict):
    for name in pyro.get_param_store():
        param_dict[name].append(pyro.param(name).item())

#Ausführung der Varianteninferenz
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

Mit ungefähr 21 Abschnitten, ungefähr 2,5 Steigungen und ungefähr 17 Dispersionen wurden vernünftige Schätzergebnisse erhalten. Aus der folgenden Abbildung ist ersichtlich, dass der Verlust und jeder Aktualisierungszielparameter ebenfalls konvergiert werden.

Beschleunigung durch GPU

Nachdem wir festgestellt haben, dass die posteriore Verteilung durch variable Inferenz geschätzt werden kann, möchten wir die Skalierbarkeit mithilfe der GPU berücksichtigen, was der größte Vorteil von Pyro ist. Der einfachste Weg, die GPU mit Pyro zu verwenden, ist am Anfang des Codes torch.set_default_tensor_type(torch.cuda.FloatTensor) Beim Schreiben werden alle generierten Tensoren von Anfang an auf die GPU gelegt. Im Beispiel betrug die Anzahl der Stichproben übrigens 100, aber wenn sich diese auf 1000, 10000, ... erhöht, wie erhöht sich dann die für die variable Inferenz erforderliche Zeit? Unter der Annahme, dass die aus 100 Stichproben geschätzten Abschnitte, Steigungen und Varianzen wahre Werte sind, werden 1000, 10000, ... Stichproben aus der wahren Verteilung (vorläufig) generiert und die Stichproben werden generiert. Wir haben die für die variable Inferenz erforderliche Zeit mit CPU und GPU als Beobachtungsdaten verglichen.

Wie Sie der Tabelle entnehmen können, wird in dem Bereich, in dem die Anzahl der Abtastwerte 10.000 oder weniger beträgt, der Vorteil der Parallelisierung nicht empfangen und die CPU ist schneller. Der Unterschied begann sich jedoch zu öffnen, als die Anzahl der Proben 100.000 überschritt, und das Ergebnis war, dass es einen Unterschied von 30 mal oder mehr für $ 10 ^ 8 $ Stücke gab. Es ist ersichtlich, dass die statistische Modellierung auch bei extrem großen Datenmengen mit GPU stressfrei durchgeführt werden kann.

Zusammenfassung

Dieses Mal stellte ich eine Methode zur Analyse durch ein einfaches Regressionsmodell unter Verwendung von Pyro vor. Der Prozess, durch den das angenommene statistische Modell Daten erzeugt, kann in der "Modell" -Methode beschrieben werden, und die posteriore Verteilung kann durch MCMC oder Differentialinferenz geschätzt werden. Wir haben auch Experimente zur Beschleunigung mit GPUs durchgeführt und bestätigt, dass Pyro skalierbare Analysen durchführen kann.