[PYTHON] Rückkehr nach dem Gaußschen Verfahren

Überblick

Hallo, das ist Kwashi. In den letzten Jahren, mit der Popularität des maschinellen Lernens, wurde es auf den industriellen Bereich angewendet, beispielsweise auf die Verwendung von Vorhersagen unter Verwendung eines Regressionsmodells zur Steuerung. Vorhersageunsicherheit ist jedoch auch in kritischen Produktkontexten wichtig. Angenommen, Sie möchten die Parameter eines Produktionsgeräts aus der Werksumgebung (Temperatur usw.) zurückgeben. Selbst wenn plötzlich ein Muster auftritt, das nicht für das Training verwendet wird, werden einige Parameter zurückgegeben. Daher ist die Unsicherheit der Inferenz auch wichtig, wie zuverlässig die Regression ist.

Die Bayes'sche Statistik ist ein Bereich, der diese Unsicherheit als Wahrscheinlichkeitsverteilung zeigt. Bayesianische Statistiken werden in verschiedenen Bereichen wie automatischem Fahren und Tonverarbeitung verwendet. Dieser Artikel enthält ein Beispiel für die Regression mithilfe der Bayes'schen Statistik. Ziehen Sie die Sündenfunktion zurück, die ein einfaches Modell als Ziel ist. Zunächst wird ein Beispiel für eine Polypoly-Regression gezeigt. Bei der Polypolyse konnte die sin-Funktion jedoch nicht ausgedrückt werden, ohne genügend Parameter einzustellen. Daher werden wir eine Methode entwickeln, um die Funktion selbst unter Verwendung des Gaußschen Prozesses abzuschätzen.

Außerdem konzentriert sich dieser Artikel auf die PyMC- und Scipy-Bibliotheken von Python.

Dieser Artikel behandelt nicht genug über Bayes'sche Statistiken. Ich habe diesen Artikel unter Bezugnahme auf die folgenden Bücher geschrieben. Bitte lesen Sie ihn. [Fortsetzung / Leicht verständliche Mustererkennung - Einführung in das Lernen ohne Lehrer](https://www.amazon.co.jp/%E7%B6%9A%E3%83%BB%E3%82%8F%E3%81%8B% E3% 82% 8A% E3% 82% 84% E3% 81% 99% E3% 81% 84% E3% 83% 91% E3% 82% BF% E3% 83% BC% E3% 83% B3% E8% AA% 8D% E8% AD% 98% E2% 80% 95% E6% 95% 99% E5% B8% AB% E3% 81% AA% E3% 81% 97% E5% AD% A6% E7% BF% 92% E5% 85% A5% E9% 96% 80% E2% 80% 95-% E7% 9F% B3% E4% BA% 95-% E5% 81% A5% E4% B8% 80% E9% 83% 8E / dp / 427421530X) Bayesian Analysis with Python

Polymere Regression der Sin-Funktion

In diesem Kapitel werden Daten gemäß einer Normalverteilung mit einer mittleren sin (x) -Standardabweichung von 0,2 unter einem Winkel x generiert, der zufällig aus einer gleichmäßigen Verteilung extrahiert wird. Es passt auch in ein Polypoly fünfter Ordnung unter Verwendung der Methode der kleinsten Quadrate. Das Polynom fünfter Ordnung ist die folgende Gleichung.

f(x) = b_0 + b_1 x + b_2 x^2 + b_3 x^3 + b_4 x^4 + b_5 x^5

Zunächst werden die Ergebnisse unten gezeigt. Die grüne Linie ist die Sinusfunktion und der blaue Punkt sind die zum Lernen verwendeten Daten. Darüber hinaus ist Rot das Ergebnis der einfachen Anpassung eines Polypolys fünfter Ordnung nach der Methode der kleinsten Quadrate.

Das Programm ist wie folgt. Die Anpassung nach der Methode der kleinsten Quadrate erfolgt mit scipy.


import numpy as np
from scipy.odr import odrpack as odr
from scipy.odr import models
import matplotlib.pyplot as plt
import theano.tensor as tt
import scipy.stats as stat
import pymc3 as pm

#Sündenfunktion
np.random.seed(5)
x = np.random.uniform(0, 10, size=20)
y = np.random.normal(np.sin(x), 0.2) #Normalverteilung, (durchschnittlich,Verteilt)

#Passform von scipy
polyFunc = models.polynomial(5) #5 ..
data = odr.Data(x, y)
myodr = odr.ODR(data, polyFunc, maxit=200)

myodr.set_job(fit_type=2) #Minimum-Quadrat-Methode

fit = myodr.run()

coeff = fit.beta[::-1]
err = fit.sd_beta[::-1]
print(coeff) #[-2.48756876e-05 -1.07573181e-02  2.09111045e-01 -1.21915241e+00, 2.11555200e+00 -1.08779827e-01] #x50

Polymerregression von Bayes

In diesem Kapitel finden wir basierend auf den im vorherigen Kapitel generierten Daten ein Polypoly fünfter Ordnung im Bayes'schen Rahmen. Das Modell entspricht dem folgenden Programm

Die t-Verteilung des Schülers

{\rm St}(x|\mu,\lambda,\nu)={\Gamma(\nu/2+1/2)\over\Gamma(\nu/2)}\left({\lambda\over\pi\nu}\right)^{1/2}\left[1+{\lambda(x-\mu)^2\over\nu}\right]^{-\nu/2-1/2}

Und. Anschließend werden die Parameter mit den Markov-Ketten-Monte-Carlo-Methoden (MCMC) gelernt.

with pm.Model() as model:
  b0 = pm.Normal('b0', mu=0, sd=10)
  b1 = pm.Normal('b1', mu=0, sd=5)
  b2 = pm.Normal('b2', mu=0, sd=5)
  b3 = pm.Normal('b3', mu=0, sd=5)
  b4 = pm.Normal('b4', mu=0, sd=5)
  b5 = pm.Normal('b5', mu=0, sd=5)

  #Cauchy-Verteilung x ∈[0, ∞) https://docs.pymc.io/api/distributions/continuous.html#pymc3.distributions.continuous.HalfCauchy
  sd = pm.HalfCauchy('sd', 5)   
  mu = b0 + b1 * x + b2 * x**2 + b3 * x**3 + b4 * x**4 + b5 * x**5
  nu = pm.Exponential('nu', 1/30)

  y_pred = pm.StudentT('y_pred', mu=mu, sd=sd, nu=nu, observed=y)
  step = pm.Metropolis()
  trace = pm.sample(100000,step=step)

Das Polypoly unter Verwendung der geschätzten Parameter ist unten rot dargestellt. Die grüne Linie ist die Sinusfunktion, der blaue Punkt sind die für das Training verwendeten Daten und der blaue Punkt ist das Ergebnis der Minimum-Square-Methode. Zusätzlich betrug der Durchschnitt der Varianzen, die als Ergebnis von MCMC erhalten wurden, 0,23, und die Varianz des erzeugten Modells konnte geschätzt werden. Es ist jedoch ersichtlich, dass der Fehler aus der tatsächlichen Sinusfunktion groß ist, und selbst in dem Teil, in dem nur wenige Daten wie etwa x = 10 vorliegen, wird die Varianz zu einem Polynom, das aus einer Konstanten 0,23 erzeugt wird. Dies ist ein schlechtes Modelldesign.

Entwicklung zum Gaußschen Prozess

Im vorherigen Kapitel wurde die Funktion f (x) durch ein polymorphes Modell dargestellt und ihre Parameter wurden abgeleitet. Die Daten sind jedoch nicht gut dargestellt und es besteht eine gute Chance, dass das Modell falsch ist. Daher verwenden wir die Methode, um direkt auf die Funktion zu schließen. Das heißt, im vorherigen Kapitel, S.(θ|x)Die Parameterschätzung wurde durchgeführt, aber in diesem Kapitel, S.(f|x)Schließen. Gaußscher Prozess als Methode, um auf diese Funktion zu schließen(GP)Benutzen.

Die offizielle Definition von GP lautet wie folgt. Alle Punkte im kontinuierlichen Raum sind Variablen zugeordnet, die einer Normalverteilung folgen, und GP ist die gleichzeitige Verteilung ihrer unendlich großen Anzahl stochastischer Variablen.

Aber ich werde hier keine besonders schwierige Diskussion führen. Es gibt drei wichtige Dinge:

  1. Die vorherige Verteilung ist eine multivariate Normalverteilung. f (x) ~ MvNormal (u | K)
  2. Die Kovarianzmatrix (K) ist der Kernel.
  3. Kovarianzmatrixparameter
  4. Die posteriore Verteilung kann analytisch berechnet werden.

Stellen Sie sich vor, der Kernel von 2 wird in einer Support-Vektor-Maschine oder dergleichen verwendet. Es gibt verschiedene Kernel, und jeder Kernel kann als Mitverteilungsmatrix verwendet werden. Dieser Artikel verwendet den Gaußschen Kernel. Werfen wir einen Blick auf die vorherige Verteilung von GP aus 1 und 2. Das folgende Programm zeichnet eine multivariate Normalverteilung (Mittelwert = 0, Kovarianzmatrix = Gaußscher Kern). Wenn Sie sich die Abbildung ansehen, sehen Sie, dass sie verschiedene Funktionen darstellt. Es gibt unendlich viele solcher Ausdrücke.

def squareDistance(x, y):
  return np.array([[(x[i]-y[j])**2 for i in range(len(x))] for j in range(len(y))])

np.random.seed(1)
tp = np.linspace(0, 10, 100)

cov = np.exp(-squareDistance(tp, tp))
plt.plot(tp, stats.multivariate_normal.rvs(cov=cov, size=6).T

Die in der obigen Abbildung gezeigte multivariate Normalverteilung wird auf verschiedene Weise durch Ändern der Kovarianzmatrix ausgedrückt. Um die Kovarianzmatrix im Bayes'schen Rahmen zu schätzen, müssen die zu schätzenden Parameter wie in 3 gezeigt eingestellt werden. Also die Kovarianzmatrix

K_{i,j} = \eta \exp(-\rho D)  \ \  (i \neq j) \\
K_{i,j} = \eta + \sigma \ \ (i = j)

Wird besorgt.

Schließlich gibt es vier analytische Berechnungen der posterioren Verteilung. Wenn die Wahrscheinlichkeit eine Normalverteilung ist, kann die GP-Posterior-Verteilung das Problem analytisch lösen, sodass die multivariate Normalverteilung, die GP-Posterior-Verteilung, mit einer einfachen Formel berechnet werden kann. Die posteriore Verteilung ermöglicht es, die Regressionsfunktion vorherzusagen.

Eine ausführliche Erläuterung der Formel finden Sie in [Gaußscher Prozess und maschinelles Lernen] Gaußscher Prozessrückgabe.

Nachfolgend finden Sie die Vorhersageergebnisse und einen Teil des Programms. Das aus der geschätzten posterioren Verteilung von GP vorhergesagte Ergebnis ist die rote Linie, die grüne Linie ist die Sinusfunktion, der blaue Punkt sind die für das Training verwendeten Daten und die blaue Linie ist das Ergebnis der Minimum-Quadrat-Methode. Aufgrund der Vorhersage durch GP liegt die rote Linie näher an der grünen Linie als die blaue Linie, sodass ersichtlich ist, dass ein Regressionsmodell mit guter Leistung direkt geschätzt werden kann, ohne die Parameter anzupassen. Wenn man sich auch ansieht, wie die roten Linien gestreut sind, kann man sehen, dass die Linien gestreut sind, das heißt, die Streuung ist groß, wenn keine Daten wie etwa 10 vorliegen.

#GP vorherige Verteilung: f(x) ~ GP(u=(u1, u2, ...), k(x, x')), u1, u2 .. =Angenommen 0
#Haftung: p(y|x, f(x)) ~ N(f, σ^2I)
#GP posteriore Verteilung p(f(x)|x,y) ~GP(U,∑)

def squareDistance(x, y):
  return np.array([[(x[i]-y[j])**2 for i in range(len(x))] for j in range(len(y))])

points = np.linspace(xmin, xmax, xsize) #Anhaltspunkt

with pm.Model() as model:
  #Vorherige Verteilung
  mu = np.zeros(xsize)
  eta = pm.HalfCauchy('eta', 5)
  rho = pm.HalfCauchy('rho', 5)
  sigma = pm.HalfCauchy('sigma', 5)

  D = squareDistance(x, x)
  K = tt.fill_diagonal(eta * pm.math.exp(-rho * D), eta + sigma) #i=eta an Position j+Speichern Sie Sigma

  gpPri = pm.MvNormal('obs', mu, cov=K, observed=y)

  #Berechnung der posterioren Verteilung(Die posteriore Verteilung der Gaußschen Verteilung kann nach der Formel berechnet werden)
  K_ss = eta * pm.math.exp(-rho * squareDistance(points, points))
  K_s = eta * pm.math.exp(-rho * squareDistance(x, points))

  MuPost = pm.Deterministic('muPost', pm.math.dot(pm.math.dot(K_s, tt.nlinalg.matrix_inverse(K)), y)) #K_s.T inv(K) y
  SigmaPost = pm.Deterministic('sigmaPost', K_ss - pm.math.dot(pm.math.dot(K_s, tt.nlinalg.matrix_inverse(K)), K_s.T))

  step = pm.Metropolis()
  trace = pm.sample(10000,step=step) 

Zusammenfassung

In diesem Artikel haben wir die Bayes'sche Statistik eingeführt, da der Ausdruck der Unsicherheit wichtig ist. Bei der Erstellung eines Regressionsmodells der Sündenfunktion war es jedoch schwierig, ein Modell zu entwerfen, das die Sündenfunktion gut ausdrückt. Daher haben wir den Gaußschen Prozess eingeführt. Ich werde detaillierte Erklärungen vermeiden und das Konzept erläutern. Weitere Informationen finden Sie in anderen Artikeln und Büchern des Artikels.

Recommended Posts

Rückkehr nach dem Gaußschen Verfahren
Gaußsche Prozessregression mit GPy
Gaußscher Prozess
PRML Kapitel 6 Gaussian Return Python-Implementierung
Gaußscher Prozess kehrt mit PyMC3 Personal Notes zurück
Gaußsche Prozessregression Numpy Implementierung und GPy
Gaußscher Prozess mit pymc3
Lineare Regressionsmethode mit Numpy
Prozess auf GPU mit chainer.cuda.elementwise
Implementieren Sie den Gaußschen Prozess in Pyro
[Statistik] [R] Versuchen Sie, die Teilungspunktregression zu verwenden.
Verarbeitung von DB-Tabelleneinfügungen mit sqlalchemy
"Gauß-Prozess und maschinelles Lernen" Gauß-Prozessregression nur mit Python-Numpy implementiert
[Für Anfänger] Prozessüberwachung mit cron
Regressionsmodell mit Scikit-Learn und dessen Visualisierung
Verwendung von Gensim mit R (Hierarchical Dirichlet Process)
Aktienkursprognose mit maschinellem Lernen (Return Edition)
[Maschinelles Lernen] Regressionsanalyse mit Scicit Learn