[PYTHON] Ein Hinweis zum Ausprobieren eines einfachen MCMC-Tutorials auf PyMC3
Überblick
Einer der kostenlosen MCMC-Sampler, die mit Python verwendet werden können, ist PyMC3. Neulich. Ein Highschool-Mädchen, das neben mir in Stava saß, sprach über Dinge wie "Vielleicht wurde es PyMC3 und schneller als PyMC2 ..." oder "Stan hat diskrete Parameter ...", also das offizielle Tutorial https: Während ich //pymc-devs.github.io/pymc3/getting_started/ ausprobierte, schrieb ich Notizen über einige Fehler und zusätzliche Experimente. Für Leute wie "Ich möchte hierarchische Felder in Python verwenden, aber es ist schwierig, MCMC mit vollem Scratch zu schreiben ..." und "Dispersive Parameter ...". Ein Student, der mit MCMC nicht sehr vertraut ist, hat es mitten in der Nacht gemacht, daher denke ich, dass es viele Teile gibt, die schwer zu lesen sind / kleine Fehler / Unverständnis, aber bitte verzeihen Sie mir.
- In MCMC ist es üblich, die Probe im frühen Stadium zu verwerfen, aber im unteren Experiment wird die Probe aufgrund der Berechnungszeit ohne Verwerfen aufgezeichnet. Daher dienen die Leistungen usw. nur als Referenz.
Umgebung
Meine PC-Umgebung ist OS X10.9.5, Python2.7 +. Mit pip können Sie eingeben, was Sie brauchen.
pip install git+https://github.com/pymc-devs/pymc3
Beachten Sie bei der Installation von PyMC3, dass das Linkziel der URL möglicherweise anders ist, wenn es sich um einen kleinen alten Artikel handelt. Es ist praktisch, Pandas und Patsy zu haben, daher denke ich, dass es besser ist, sie mit pip einzulegen.
pip install pandas
pip install patsy
Grundlegende Verarbeitung
Da ich einen Teil des Tutorials verwenden werde, um zu erklären, wo der Fehler aufgetreten ist, werde ich die minimal erforderlichen Teile aus der ersten Hälfte des Tutorials zusammenfassen. In diesem Artikel wird der PyMC3-abhängige Teil explizit in Form von pm geschrieben. Die Uhr wird im Tutorial weggelassen, aber der Code ist bis auf diesen derselbe.
Zunächst etwas Import.
import numpy as np
%pylab inline
np.random.seed(27)
import pymc3 as pm
Erstellen Sie einen Beispieldatensatz wie folgt.
alpha, sigma = 1, 1
beta = [1, 2.5]
size = 100
X1 = np.linspace(0, 1, size)
X2 = np.linspace(0,.2, size)
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma
Definieren Sie das Modell.
#Ein Modell, das Y auf X1 und X2 linear zurückgibt
#Bedeutet, ein Modell mit dem Namen model zu definieren
with pm.Model() as model:
# Priors for unknown model parameters
alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10, shape=2)
sigma = pm.HalfNormal('sigma', sd=1)
# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2
# Likelihood (sampling distribution) of observations
Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
Führen Sie MCMC mit dem obigen Modell aus. (Mit meiner CPU hat es 13,6 Sekunden gedauert)
from scipy import optimize
with model:
# obtain starting values via MAP
start = pm.find_MAP(fmin=optimize.fmin_powell)
# instantiate sampler
step = pm.Slice(vars=[sigma])
# draw 5000 posterior samples
trace = pm.sample(5000, start=start, step=step)
Für das Argument von pm.sample () ist die erste 5000 die Größe der Stichprobe, die Sie erhalten möchten. Beim nächsten Start wird der Anfangswert gesetzt. Dieses Mal wird die MAP-Schätzung verwendet, es ist jedoch möglich, sie nicht anzugeben. Die Schnittabtastung wird im letzten Schritt angegeben, aber wenn sie nicht angegeben wird, wird sie angegeben. ・ Binäre Variable: Binäre Metropole ・ Diskrete Variable: Metropolis ・ Kontinuierliche Variable: NUTS Ist standardmäßig zugewiesen. In dem hier verwendeten einfachen Modell gab es je nach Stichprobenmethode fast keinen Unterschied. Der Hamiltonian MC und NUTS, die mit Stan verwendet werden können, können auch mit PyMC3 verwendet werden.
Die Visualisierung ist sehr einfach, nur eine Aufnahme mit Traceplot wie folgt.
pm.traceplot(trace[4000:])
Die zusammenfassenden Informationen können auch wie folgt angezeigt werden.
pm.summary(trace[4000:])
Der Einfachheit halber wurde diesmal nur die Normalverteilung behandelt, aber natürlich können verschiedene andere Verteilungen wie Gammaverteilung, Betaverteilung, Binomialverteilung, Poissonverteilung usw. behandelt werden.
- Für die Beziehung zu PyMC2 war http://breakbee.hatenablog.jp/entry/2014/08/04/031342 usw. hilfreich.
Ärger
Obwohl es sich um ein so praktisches PyMC3 handelt, gibt es einige Probleme. Versuchen wir das Tutorial.
1. Ich erhalte eine Fehlermeldung, wenn ich die Methode theano verwende
Ich glaube, ich habe neulich gerade eine Diskussion darüber auf Twitter gesehen, aber in diesem Fall ist die Fehlerursache wahrscheinlich die letzte von https://github.com/pymc-devs/pymc3/blob/master/pymc3/model.py Linie
theano.config.compute_test_value = 'raise'
(Referenz: https://sites.google.com/site/iwanamidatascience/vol1/support_tokushu). Ich kenne die Details des Caches nicht, aber nach dem Import von PyMC3
import theano
theano.config.compute_test_value = 'ignore'
Ich konnte es lösen, indem ich mein Herz eintippte. Übrigens scheint es, dass Sie anstelle von "Signore" "abschreiben" können, aber in glm Berechnung usw.
scratchpad instance has no attribute 'test_value'
Ich denke, es ist sicher, es vorerst auf "Signore" zu setzen. Übrigens ist es auch auf der Seite Probabilistic Matrix Factorization der offiziellen Beispiele zu finden.
2. Problematischer Umgang mit selbst erstellten Funktionen und Klassen
2.1 as_op unnötige Theorie
Im Abschnitt Arbitrary Deterministics des Tutorials wird die deterministische Variable b durch die folgende Funktion crazy_modulo3 () gefunden.
import theano.tensor as T
from theano.compile.ops import as_op
@as_op(itypes=[T.lscalar], otypes=[T.lscalar])
def crazy_modulo3(value):
if value > 0:
return value % 3
else :
return (-value + 1) % 3
with pm.Model() as model_deterministic:
a = pm.Poisson('a', 1)
b = crazy_modulo3(a)
Der Teil @ as_op () ist wie ein Zauberspruch für die Verwendung einer Funktion, die eine Variable vom Typ theano.tensor berechnet. Er wird unmittelbar vor der Definition der Funktion zum Deklarieren des Eingabe- / Ausgabetyps in die Zeile geschrieben. Dies wird auch im PyMC3-Tutorial beschrieben.
Theano needs to know the types of the inputs and outputs of a function, which are specified for as_op by itypes for inputs and otypes for outputs.
** Dies kann jedoch tatsächlich wie folgt definiert werden, ohne @ as_op zu verwenden. ** **.
import theano.ifelse
def symbolf(x):
return theano.ifelse.ifelse(T.gt(x,0), x%3, (-x+1)%3)
def crazy_modulo3(value):
return symbolf(value)
Ich bin ein bisschen süchtig danach, also werde ich es erklären.
Zum Beispiel naiv ohne ifelse
def crazy_modulo3(value):
if value > 0:
return value % 3
else :
return (-value + 1) % 3
Dann werde ich wütend, wenn ich "Wert> 0" nicht vergleichen kann. Ebenfalls,
def crazy_modulo3(value):
if T.lt(0,value):
return value % 3
else :
return (-value + 1) % 3
Dann tritt kein Fehler auf, aber wenn ich den Wert verfolge, wird "T.lt (0, Wert)" nicht falsch (Theano hat keinen booleschen Typ. Stattdessen werden alle booleschen Tensoren in "int8" dargestellt. . --http: //deeplearning.net/software/theano/library/tensor/basic.html).
Daher wird ifelse zum erzwungenen Vergleich verwendet. Dies ist ein bisschen schmerzhaft, aber jetzt kann ich die Funktion definieren, ohne @ as_op zu verwenden.
2.2 Verfolgen von deterministischen Plotvariablen
Obwohl es in Ordnung ist, die deterministische Variable b im obigen Code zu definieren, beschreibt das Lernprogramm die nachfolgende Stichprobenmethode nicht. Probieren Sie es dort aus
with model_deterministic:
trace = pm.sample(500)
pm.traceplot(trace)
Wenn Sie es versuchen, wird nur der Wert der Variablen a dargestellt. Ich war auch an der Variablen b interessiert, also fand ich beim Betrachten des Inhalts von pm. So etwas wie "pm.Deterministic ()". Der Rest ist die Definition des Modells
b = crazy_modulo3(a)
Zu
b = pm.Deterministic("b",crazy_modulo3(a))
Wenn Sie es ändern und abtasten, sollte der Wert von b korrekt dargestellt werden.
2.3 Vermeiden Sie die Verwendung von as_op zur Berechnung von Zufallsvariablen
Im Abschnitt Beliebige Verteilungen des Lernprogramms wird wie folgt beschrieben, wie Sie die Wahrscheinlichkeitsverteilung definieren, die Sie abtasten möchten.
with pm.Model() as model1:
alpha = pm.Uniform('intercept', -100, 100)
# Create custom densities
beta = pm.DensityDist('beta', lambda value: -1.5 * T.log(1 + value**2), testval=0)
eps = pm.DensityDist('eps', lambda value: -T.log(T.abs_(value)), testval=1)
# Create likelihood
like = pm.Normal('y_est', mu=alpha + beta * X, sd=eps, observed=Y)
Diese durch den Lambda-Ausdruck und "pm.DensityDist" definierte Methode ist kein Problem. Wenn Sie die Probenahme wie folgt fortsetzen, funktioniert sie ordnungsgemäß.
with model1:
trace = pm.sample(500)
pm.traceplot(trace)
Als Code mit ähnlicher Bedeutung verwendet das Lernprogramm die folgende Methode zum Definieren einer Klasse (** Ich werde darauf hinweisen, dass dieser Code wahrscheinlich bald falsch sein wird **). Dies verwendet keine Lambda-Ausdrücke, daher scheint es den Vorteil zu haben, das Schreiben von Funktionen zu vereinfachen. Mit Ausnahme von Beta sind dies die Parameter, die ich hinzugefügt habe, damit sie die gleiche Form wie Modell 1 oben haben.
class Beta(pm.distributions.Continuous):
def __init__(self, mu, *args, **kwargs):
super(Beta, self).__init__(*args, **kwargs)
self.mu = mu
self.mode = mu
def logp(self, value):
mu = self.mu
return beta_logp(value - mu)
@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
return -1.5 * np.log(1 + (value)**2)
with pm.Model() as model2:
beta = Beta('slope', mu=0, testval=0)
#I add other parameters to follow model1 above
alpha = pm.Uniform('intercept', -100, 100)
eps = pm.DensityDist('eps', lambda value: -T.log(T.abs_(value)), testval=1)
like = pm.Normal('y_est', mu=alpha + beta * X, sd=eps, observed=Y)
Das Beispiel "@ as_op" erscheint in der Beta-Berechnung. Als ich damit früher die deterministische Variable berechnet habe, hat die Stichprobe gut funktioniert, aber was ist mit dieser Zeit?
with model2:
trace = pm.sample(100)
pm.traceplot(trace)
Wenn Sie dies tun, erhalten Sie wahrscheinlich einen ↓ Fehler.
AttributeError: 'FromFunctionOp' object has no attribute 'grad'
Anscheinend hat die von ** @ as_op eingebrachte Nachricht nicht das Gradientenberechnungsattribut grad, das für die MCMC-Abtastung ** erforderlich ist. War es in Ordnung, weil der Gradient nicht mit deterministischen Variablen berechnet wird? Als ich hier googelte, fand ich eine Antwort von einer großartigen Person: "Ich weiß nicht, wie ich damit umgehen soll."
The only way to get the automatic gradient computation is by expressing your density in terms of theano operators. as_op creates a blackbox function for which autodiff will not work so there is no way I know of (except numerical differentiation) to make this work. https://github.com/pymc-devs/pymc3/issues/601
Hier ist jedoch das Ergebnis von 2.1. Mit anderen Worten, wenn @ as_op schlecht ist, warum nicht zuerst ohne @ as_op definieren? Es ist eine Geschichte. Deshalb,
@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
return -1.5 * np.log(1 + (value)**2)
Wird wie folgt umgeschrieben.
def beta_logp(value):
return -1.5 * T.log(1 + (value)**2)
Wie in 2.1 können Sie sicher rechnen, wenn Sie beim Umgang mit Variablen vom Typ theano.tensor die Regeln beachten (in diesem Fall verwenden Sie T.log anstelle von np.log). Das Ergebnis ist fast das gleiche (obwohl es einige Unschärfen gibt, weil es sich um Stichproben handelt).
** Wenn Sie Ihre eigenen Funktionen und Klassen verwenden, unabhängig davon, ob Sie Lambda-Ausdrücke verwenden oder nicht, müssen Sie beim Codieren auf die Syntax von theano achten. ** **.
3. Geschwindigkeit und Leistung bei glm
Abschließend werde ich auf die Geschwindigkeit und Leistung bei der Berechnung des verallgemeinerten linearen Modells (glm) eingehen. Details finden Sie unter https://github.com/pymc-devs/pymc3/issues/544, aber "NUTS ist oft standardmäßig langsam" "Metropolis ist schnell" "Abhängig von der Art der Daten und der Maschine, wenn Hamiltonian MC schnell ist" Es gibt auch eine Geschichte. " Als Referenz habe ich die Tutorial-Daten auf glm (gewöhnliche lineare Regression) in meiner Umgebung angewendet.
import pandas as pd
df = pd.DataFrame({'x1': X1, 'x2': X2, 'y': Y})
with pm.Model() as model:
pm.glm.glm('y ~ x1+x2', df)
start = pm.find_MAP(fmin=optimize.fmin_powell)
step = pm.HamiltonianMC()
trace = pm.sample(1000, start=start, step=step)
pm.traceplot(trace)
・ HamiltonianMC: In 1,7 Sekunden fertig. Konvergenz ist auch in Ordnung. (Abbildung 1) ・ Metropole: Fertig in 0,3 Sekunden. Konvergenz ist nicht gut. ・ NÜSSE: Selbst nach 5 Minuten können zu spät nur 8 Proben entnommen werden -Handgeschriebene Probenahme ohne Verwendung der in 1. dieses Artikels durchgeführten pm.glm-Klasse: Fertig in 3,2 Sekunden. Konvergenz ist Hamiltonian MC unterlegen, aber nicht so schlecht. (Figur 2)
Wahre Parameter: "Alpha = 1, Beta [0] = 1, Beta [1] = 2,5, Sigma = 1"
(Abbildung 1)
(Abb. 2: Die blaue Linie von Beta entspricht dem Regressionskoeffizienten von X1, die grüne Linie entspricht dem Regressionskoeffizienten von X2 und Sigma entspricht sd)
Das Ergebnis war. Übrigens, wenn ich NUTS in der Mitte abschneide und es zeichne, scheint die Skalierung der Parameter seltsam zu sein. Als Lösung in diesem Fall auf der obigen Problemseite
C = pm.approx_hessian(model.test_point)
step = pm.NUTS(scaling=C)
Wenn ja, steht geschrieben, dass die Geschwindigkeit herauskommt. Als ich die Version von Sphinx aktualisierte (wenn Sie sie nicht haben, installieren Sie sie bitte mit pip) und die obigen zwei Zeilen eintippte, wurde mir gesagt, dass das führende Moll kein positiver Wert ist und nicht berechnet werden kann.
LinAlgError: 2-th leading minor not positive definite
Dies scheint ein wenig unlösbar zu sein, also habe ich es durchgearbeitet, aber hier scheint es eine verwandte Diskussion zu geben → https://github.com/scikit-learn/scikit-learn/issues/2640
Als persönlicher Eindruck: ** Wenn Sie nur glm verwenden möchten, versuchen Sie es mit glm by NUTS und verwenden Sie approx_hessian für die Skalierung. Wenn dies nicht funktioniert, versuchen Sie es mit glm by HMC. ** **.
Außerdem scheint die Methode, die die pm.glm-Klasse nicht verwendet, wie im Beispiel von 1 beschrieben, die Möglichkeit zu haben, die Leistung zu verbessern, wenn Sie beispielsweise den vorherigen Brunnen festlegen und das Verhalten des Hyperparameters oder sehen möchten Wenn ich ein komplizierteres hierarchisches Bayes'sches Modell verwenden möchte, dachte ich, ich sollte dies übernehmen, anstatt pm.glm zu verwenden.
Schließlich
Ich habe es lange Zeit mit meinem eigenen Memo geschrieben, aber ich hoffe, es wird eine kleine Referenz für den Ort sein, der süchtig zu machen scheint. Bei @ as_op gab es etwas, das ziemlich umständlich war, wie zum Beispiel eine Nacht damit zu verschwenden, den Grad auf der Theano-Seite neu zu definieren.
Außerdem habe ich auf der offiziellen Seite das Gefühl, dass Beispiele interessanter sind als das Tutorial (Probabilistische Matrixfaktorisierung, Überlebensanalyse usw.), daher werde ich es erneut versuchen, wenn ich mich wie Mitternacht fühle.
- Die Definition von grad finden Sie unter http://deeplearning.net/software/theano/extending/op.html. Die Eingabe, wenn Sie den ursprünglich definierten T.grad als eine Funktion betrachten, die mit theano.function erstellt werden kann, scheint eine Liste zu sein. Wenn Sie jedoch diesem Lernprogramm folgen, können Sie keine Funktion definieren, deren Eingabe eine Liste ist, und aufgeben tat.
Ich werde es entsprechend hinzufügen.
Recommended Posts