[Python] Bayesianische Schätzung mit Pyro

Ich habe kürzlich von Pyro erfahren, das mit Wahrscheinlichkeitsmodellen umgeht, und dachte, es wäre interessant, also habe ich es als Versuch versucht. Dieser Artikel teilt den Quellcode. Der Quellcode ist in Jupyter Notebook geschrieben. Bitte beachten Sie, dass es fast keine theoretische Erklärung gibt.

Wir veröffentlichen auch die ipynb-Datei auf github. https://github.com/isuya0508/exercise_Pyro/blob/master/bayesian_inference_bernulli.ipynb

Was ist Pyro?

Pyro ist eine Bibliothek, die probabilistische Modelle mit Pytorch als Backend verarbeitet. Sie können es von pip installieren.

pip install pyro-ppl 

Zu diesem Zeitpunkt muss Pytorch im Voraus installiert werden. Weitere Informationen finden Sie auf der offiziellen Seite.

Versuch es

Lassen Sie uns diesmal den Parameter $ p $ aus den Daten schätzen, die der Bernoulli-Verteilung $ Ber (p) $ folgen. Importieren Sie zunächst die erforderlichen Module.

from collections import Counter

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

import torch
import torch.distributions.constraints as constraints
import pyro
import pyro.distributions as dist
from pyro.optim import SGD, Adam
from pyro.infer import SVI, Trace_ELBO

%matplotlib inline

pyro.set_rng_seed(0)
pyro.enable_validation(True)

Datengenerierung

Erstellen Sie Daten mit Zufallszahlen. Zu diesem Zeitpunkt muss der Typ der Pytorch-Tensor sein.

obs = stats.bernoulli.rvs(0.7, size=30, random_state=1)
obs = torch.tensor(obs, dtype=torch.float32)
obs

> tensor([1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1.,
          1., 1., 0., 0., 1., 1., 0., 0., 1., 1., 1., 0.])

Überprüfen Sie die Anzahl der Einsen in den Daten.

Counter(obs.numpy())

> Counter({1.0: 23, 0.0: 7})

Daher ist die wahrscheinlichste Schätzung für den Parameter $ p $ $ 23/30 \ fallenddotseq 0,77 $. [^ 2] Danach werden wir Bayesian $ p $ schätzen.

Festlegen der vorherigen Verteilung

Bei der Bayes'schen Schätzung wird die vorherige Verteilung der Parameter angenommen und die hintere Verteilung wird durch Kombinieren der beobachteten Daten erhalten.

Für den Parameter $ p $ der Bernoulli-Verteilung $ Ber (p) $ wird häufig eine Beta-Verteilung als vorherige Verteilung angenommen [^ 1] In Pyro beschreibt die "Modell" -Methode die vorherige Verteilung und das Modell der Daten.

def model(data):
    #Vorausgesetzte Verteilung voraussetzen
    prior = dist.Beta(1, 1)

    #Datenmodellierung
    p = pyro.sample('p', prior)
    for i in range(len(data)):
        pyro.sample(f'obs{i}', dist.Bernoulli(p), obs=data[i])

Diesmal nehmen wir $ Beta (1, 1) $ an, was tatsächlich mit einer gleichmäßigen Verteilung übereinstimmt.

Hintere Verteilung einstellen

Beschreiben Sie die posteriore Verteilung mit der "Guide" -Methode. Die hintere Verteilung sollte sowohl eine Beta-Verteilung als auch die vorherige Verteilung sein. Geben Sie zu diesem Zeitpunkt einen geeigneten Anfangswert als Parameter für die posteriore Verteilung an.

def guide(data):
    #Definition der posterioren Verteilung
    alpha_q = pyro.param('alpha_q', torch.tensor(15), constraint=constraints.positive)
    beta_q = pyro.param('beta_q', torch.tensor(15), constraint=constraints.positive)
    posterior = dist.Beta(alpha_q, beta_q)

    pyro.sample('p', posterior)

Dieser hintere Verteilungsparameter $ \ alpha, \ beta $ wird gefunden.

Nachverteilte Anpassung

Für die Methode zur Schätzung der Parameter der posterioren Verteilung verwenden wir diesmal die probabilistische Variationsschätzung. Es scheint, dass Pyro diese Methode grundsätzlich verwendet. In diesem Beispiel der Bernoulli-Verteilung ist es Unsinn, eine Approximationsmethode wie die Variationsschätzung zu verwenden, da die posteriore Verteilung analytisch erhalten werden kann, aber ich werde diese Methode für die Praxis verwenden. (Ursprünglich ist dies die Methode, die verwendet wird, wenn die Verteilung nicht analytisch berechnet werden kann.)

Generieren Sie für die Implementierung eine Instanz von Gradient Descent (SGD) und eine Instanz von Probabilistic Variant Estimation (SVI). Wenn Sie eine SVI-Instanz erstellen, geben Sie ihr das oben definierte "Modell" und "Handbuch".

optimizer = SGD(dict(lr=0.0001, momentum=0.9))
svi = SVI(model, guide, optimizer, loss=Trace_ELBO())

Jetzt müssen Sie nur noch "svi.step (obs)" wiederholen, um die posterioren Verteilungsparameter zu aktualisieren.

NUM_STEPS = 2000
pyro.clear_param_store()

history = {
    'loss': [],
    'alpha': [],
    'beta': []
}

for step in range(1, NUM_STEPS + 1):
    loss = svi.step(obs)
    
    history['loss'].append(loss)
    history['alpha'].append(pyro.param('alpha_q').item())
    history['beta'].append(pyro.param('beta_q').item())
    
    if step % 100 == 0:
        print(f'STEP: {step} LOSS: {loss}')

> 
STEP: 100 LOSS: 17.461310371756554
STEP: 200 LOSS: 18.102468490600586
(Unterlassung)
STEP: 1900 LOSS: 17.97793820500374
STEP: 2000 LOSS: 17.95139753818512

Hier werden der Verlust, der angepasst wird, und die posterioren Verteilungsparameter $ \ alpha und \ beta $ in der "Geschichte" aufgezeichnet. Die Darstellung dieser Zahlen für jeden Schritt sieht folgendermaßen aus: (Der Quellcode wird weggelassen.)

!

Überprüfen Sie das endgültig erhaltene $ \ alpha, \ beta $ und den erwarteten Wert und die Varianz der posterioren Verteilung.

infered_alpha = pyro.param('alpha_q').item()
infered_beta = pyro.param('beta_q').item()
posterior = stats.beta(infered_alpha_beta, infered_beta_beta)

print(f'alpha: {infered_alpha}')
print(f'beta: {infered_beta}')
print(f'Expected: {posterior.expect()}')
print(f'Variance: {posterior.var()}')

>
alpha: 25.764650344848633
beta: 7.556574821472168
Expected: 0.7732203787899605
Variance: 0.005109101547631603

Zeichnen wir die Vor- und Nachverteilung. (Der Quellcode wird weggelassen.)

Es scheint, dass es gut geschätzt werden kann.

Zusammenfassung

Ich habe versucht, mit Pyro eine einfache Bayes'sche Schätzung durchzuführen. Diesmal war es eine einfache Bayes'sche Schätzung, aber ich denke, Pyros Stärke ist, dass es komplexe Modelle flexibel und einfach beschreiben kann. Es gibt viele Beispiele für probabilistische Methoden in der offiziellen Dokumentation, und es kann eine Lernerfahrung sein, sie nur anzusehen.

[^ 2]: Die wahrscheinlichste Schätzung von $ p $ entspricht dem Stichprobenmittelwert. [^ 1]: Es ist bekannt, dass die hintere Verteilung auch der Beta-Verteilung folgt, indem die vorherige Verteilung von $ p $ zu einer Beta-Verteilung gemacht wird. Eine solche Verteilung wird als konjugierte vorherige Verteilung bezeichnet.