[PYTHON] Gaußsche Prozessregression Numpy Implementierung und GPy

Einführung

Möglicherweise sehen Sie ein Diagramm, in dem Sie sich fragen, ob Sie eine solche Regressionslinie zeichnen können. Zum Beispiel

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Times New Roman"      #Stellen Sie die gesamte Schriftart ein
plt.rcParams["xtick.direction"] = "in"               #Nach innen die x-Achsen-Skalierungslinie
plt.rcParams["ytick.direction"] = "in"               #Y-Achsen-Skalierungslinie nach innen
plt.rcParams["xtick.minor.visible"] = True           #Hinzufügen einer x-Achsen-Hilfsskala
plt.rcParams["ytick.minor.visible"] = True           #Hinzufügen einer Hilfsskala der y-Achse
plt.rcParams["xtick.major.width"] = 1.5              #Linienbreite der x-Achsen-Hauptskalenlinie
plt.rcParams["ytick.major.width"] = 1.5              #Linienbreite der Hauptskalenlinie der y-Achse
plt.rcParams["xtick.minor.width"] = 1.0              #Linienbreite der Hilfsskalenlinie der x-Achse
plt.rcParams["ytick.minor.width"] = 1.0              #Linienbreite der Hilfsskalenlinie der y-Achse
plt.rcParams["xtick.major.size"] = 10                #Länge der x-Achsen-Hauptskalenlinie
plt.rcParams["ytick.major.size"] = 10                #Länge der Hauptmaßstablinie der y-Achse
plt.rcParams["xtick.minor.size"] = 5                 #Länge der Hilfsskala der x-Achse
plt.rcParams["ytick.minor.size"] = 5                 #Länge der Hilfsskala der y-Achse
plt.rcParams["font.size"] = 14                       #Schriftgröße
plt.rcParams["axes.linewidth"] = 1.5                 #Gehäusedicke

def func(x):
    return 1.5*x+2.0 + np.random.normal(0.0,0.1,len(x))

x = np.array([0.0,0.2,0.5,1.0,10.0])
y = func(x) #Messpunkt
a,b = np.polyfit(x,y,1) #Lineare Regression

fig,axes = plt.subplots()
axes.scatter(x,y)
axes.plot(x,a*x+b)

Durch lineare Regression konnte ich eine gerade Linie zeichnen. Es gibt jedoch keine Messpunkte zwischen 1 und 10, und dieser Teil scheint intuitiv unsicher zu sein. Als ich nach einer guten Regressionsanalyse suchte, die eine solche Unsicherheit ausdrücken könnte, fand ich eine sogenannte Gaußsche Prozessregression, daher möchte ich sie für Studien implementieren und ausprobieren.

Referenz

Sogar ich, der ich mit Mathematik nicht vertraut bin, konnte es lesen. Dieser Beitrag entspricht den Kapiteln 2 und 3. Gaußscher Prozess und maschinelles Lernen

Hyperparameter existierten in der Gaußschen Prozessregression, und die MCMC-Methode wurde verwendet, um sie abzuschätzen. Probenahme nach der Markov-Ketten-Monte-Carlo-Methode (MCMC) mit dem Moderator

Es gibt eine Bibliothek, ohne sie selbst implementieren zu müssen. [Gpy vs scikit-learn: Gaußscher Prozess kehrt mit Python zurück](https://nykergoto.hatenablog.jp/entry/2017/05/29/python%E3%81%A7%E3%82%AC%E3%82% A6% E3% 82% B9% E9% 81% 8E% E7% A8% 8B% E5% 9B% 9E% E5% B8% B0_ ~ % E3% 83% A2% E3% 82% B8% E3% 83% A5% E3% 83% BC% E3% 83% AB% E3% 81% AE% E6% AF% 94% E8% BC% 83 ~)

Gaußsche Prozessregression

1. Gaußscher Prozess

Für jede natürliche Zahl $ N $ entspricht der Ausgabevektor der Eingabe $ x_1, x_2, ..., x_N $ ${\mathbf f}=(f(x_1),f(x_2),...,f(x_N))$ Ist der Durchschnitt $ {\ mathbf \ mu} = (\ mu (x_1), \ mu (x_2), ..., \ mu (x_N)) $, $ K_ {nn ^ {'}} = k (x_n, x_n) ^ {'}) Wenn Sie der Gaußschen Verteilung $ {\ mathbf \ mu} folgen, {\ rm K}) $ mit der Matrix $ \ rm K $ als Element und $ {\ mathcal N} ({\ mathbf \ mu}, {\ rm K}) $ als Kovarianzmatrix, $ f $ Soll dem Gaußschen Prozess folgen $f \sim {\rm GP}(\mu({\mathbf x}), k({\mathbf x},{\mathbf x^{'}}))$ Schreiben.

Ich bin mir jedoch nicht sicher über die Art des Gaußschen Prozesses, daher werde ich versuchen, die Kernelmatrix unter Verwendung der radialen Basisfunktion als Kernelfunktion und Stichprobe aus der Gaußschen Verteilung zu berechnen, die der Kernelmatrix folgt.

def rbf(x,s=[1.0,1.0,1.0e-6]):
    """
    Radial Basis Function, RBF :Radiale Basisfunktion
    """
    s1,s2,s3 = s[0],s[1],s[2]
    
    X, X_ = np.meshgrid(x,x)
    K = s1**2 * np.exp(-((X-X_)/s2)**2) + np.identity(len(x)) * s3**2
    return K

N = 100
x = np.linspace(0.0,10.0,N) #Eingang
K = rbf(x) #Kernel-Matrix

fig,axes = plt.subplots()
axes.imshow(K)

Um aus einer mehrdimensionalen Gaußschen Verteilung mit einer verteilten, gemeinsam verteilten Matrix von $ \ rm K $ eine Stichprobe zu erstellen, generieren Sie zufällig $ \ mathbf x $ und konvertieren Sie sie in $ \ mathbf {y} = \ mathbf {Lx} $ .. $ \ mathbf L $ kann durch Colleskey-Zerlegung von $ {\ rm K} $ erhalten werden.

L = np.linalg.cholesky(K) #Colleskey-Zerlegung der Kernelmatrix

fig,axes = plt.subplots()
axes.plot()
for i in range(10):
    _x = np.random.normal(0.0,1.0,N)
    axes.plot(x,np.dot(L,_x))

Funktionen werden zufällig angezeigt, genau wie Zahlen, wenn Sie würfeln. Es scheint auch, dass die intuitive Natur des Gaußschen Prozesses darin besteht, dass ähnliche Werte aus ähnlichen Werten hervorgehen.

2. Gaußsche Prozessregression

Es gibt $ N $ Messungen. Der Einfachheit halber ist der Durchschnittswert von y 0. ${\mathcal D} = \{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\}$ Zu diesem Zeitpunkt zwischen $ \ mathbf x $ und $ y $ $y=f({\mathbf x})$ Diese Funktion $ f $ ist ein Gaußscher Prozess mit einem Durchschnitt von 0 $f \sim {\rm GP}(0,k(x,x^{'}))$ Angenommen, es wird generiert aus. Wenn Sie $ {\ mathbf y} = (y_1, y_2, ..., y_N) ^ {T} $ sagen, folgt dieses $ \ mathbf y $ einer Gaußschen Verteilung und der Kernelfunktion $ k $ für jedes Eingabepaar. Gebrauch machen von ${\rm K}\ =k(x,x^{'})$ Verwenden der in angegebenen Kernelmatrix $ \ rm K $ $y \sim {\mathcal N}(0,{\rm K})$ Ist festgelegt.

So finden Sie den Wert von $ y ^ {\ ast} $ bei $ x ^ {\ ast} $, der nicht in den Daten enthalten ist. $ \ Mathbf y $ plus $ y ^ {\ ast} $ $ {\ mathbf y} ^ {'} = (y_1, y_2, ..., y_N, y ^ {\ ast}) ^ { Sei T} $. Wenn Sie $ \ rm K ^ {'} $ in der aus $ (x_1, x_2, ..., x_N, x ^ {\ ast}) ^ {T} $ berechneten Kernelmatrix ausführen, folgen alle ebenfalls der Gaußschen Verteilung. Damit $y^{'} \sim {\mathcal N}(0,{\rm K}^{'})$ Es wird sein. Das ist, $\begin{pmatrix} {\mathbf y} \\\ y^{\ast} \end{pmatrix} \sim {\mathcal N} \begin{pmatrix}0, {\begin{pmatrix}{\rm K} & {k_{\ast}} \\\ {k_{\ast}^T} & {k_{\ast \ast}}\end{pmatrix}}\end{pmatrix}$ Ist festgelegt. Da dies ein Ausdruck für die gleichzeitige Verteilung von $ y ^ {\ ast} $ und $ \ mathbf y $ ist, ist die bedingte Wahrscheinlichkeit von $ y ^ {\ ast} $ bei $ \ mathbf y $ eine Gaußsche Verteilung. Sie wird aus der bedingten Wahrscheinlichkeit zwischen den Elementen von berechnet. Die bedingte Wahrscheinlichkeit ist: $p(y^{\ast}|x^{\ast},{\mathcal D}) = {\mathcal N}(k_{\ast}^T {\rm K}^{-1} {\mathbf y}, k_{\ast \ast}-k_{\ast}^T {\rm K}^{-1} k_{\ast})$

Schreiben Sie den Code gemäß der obigen Formel.

def rbf(x_train,s,x_pred=None):
    """
    Radial Basis Function, RBF :Radiale Basisfunktion
    """
    s1,s2,s3 = s[0],s[1],s[2]

    if x_pred is None:
        x = x_train
        X, X_ = np.meshgrid(x,x)
        K = s1**2 * np.exp(-((X-X_)/s2)**2) + np.identity(len(x)) * s3**2 #Kernel-Matrix
        return K
    else:
        x = np.append(x_train,x_pred)
        X, X_ = np.meshgrid(x,x)
        K_all = s1**2 * np.exp(-((X-X_)/s2)**2) + np.identity(len(x)) * s3**2 #Kernel-Matrix
        K = K_all[:len(x_train),:len(x_train)]
        K_s = K_all[:len(x_train),len(x_train):]
        K_ss = K_all[len(x_train):,len(x_train):]
        return K,K_s,K_ss
    
def pred(x_train,y_train,x_pred,s):
    K,K_s,K_ss = rbf(x_train,s,x_pred=x_pred)
    K_inv = np.linalg.inv(K) #Inverse Matrix
    y_pred_mean = np.dot(np.dot(K_s.T,K_inv), y_train) #Erwarteter Wert von y
    y_pred_cov = K_ss - np.dot(np.dot(K_s.T,K_inv), K_s) #Verteilte mitverteilte Matrix
    y_pred_std = np.sqrt(np.diag(y_pred_cov)) #Standardabweichung
    return y_pred_mean,y_pred_std

Erstellen Sie eine geeignete Funktion und probieren Sie sie als Messwert aus.

def func(x):
    return np.sin(2.0*np.pi*0.2*x) + np.sin(2.0*np.pi*0.1*x)

pred_N = 100 #Anzahl der vorhergesagten Punkte
N = 10 #Anzahl der Messpunkte
x_train = np.random.uniform(0.0,10.0,N) #Messpunkt:x
y_train = func(x_train) + np.random.normal(0,0.1,1) #Messpunkt:y

x_true = x_pred = np.linspace(-2.0,12.0,pred_N)
fig,axes = plt.subplots()
axes.plot(x_true, func(x_true), label="True")
axes.scatter(x_train,y_train, label="Measured")
axes.legend(loc="best")

Lassen Sie uns eine Gaußsche Prozessregression durchführen.

s = [1.0,1.0,0.1]
y_pred_mean,y_pred_std = pred(x_train,y_train,x_pred,s)

fig,axes = plt.subplots(figsize=(8,6))
axes.set_title("s1=1.0, s2=1.0, s3=0.1")
axes.plot(x_true, func(x_true), label="True")
axes.scatter(x_train,y_train, label="Measured")
axes.plot(x_pred, y_pred_mean, label="Predict")
plt.fill_between(x_pred,y_pred_mean+y_pred_std,y_pred_mean-y_pred_std,facecolor="b",alpha=0.3)
axes.legend(loc="best")
axes.set_xlim(-2.0,12.0)

Ich konnte so zurückkehren. Die Verteilung ist weit verbreitet, wenn die Anzahl der Messpunkte gering ist. Nun werden die Hyperparameter der Kernelfunktion von Hand angegeben. Mal sehen, was passiert, wenn wir das ändern.

Das Ergebnis hat sich sehr verändert. Ich möchte auch die Hyperparameter schätzen.

Setzen Sie die Hyperparameter als $ \ theta $ zusammen. Die Kernelmatrix hängt von $ \ theta $ ab, wo die Wahrscheinlichkeit von Trainingsdaten liegt

Wenn Sie das Protokoll nehmen $\ln{p({\mathbf y}|{\mathbf x},\theta)} \propto -\ln{|{\rm K}|}-{\mathbf y}^T {\rm K}_{\theta}^{-1} {\mathbf y} + const.$ Und findet $ \ theta $, das die obige Gleichung maximiert.

Die MCMC-Methode wird verwendet, da es bei der Gradientenmethode leicht ist, in der lokalen Lösung hängen zu bleiben.

import emcee
def objective(s): #Zielfunktion
    K = rbf(x_train,s)
    return -np.linalg.slogdet(K)[1] - y_train.T.dot(np.linalg.inv(K)).dot(y_train) #Beweise

ndim = 3
nwalker = 6
s0 = np.random.uniform(0.0,5.0,[nwalker,ndim]) #Ausgangsposition
sampler = emcee.EnsembleSampler(nwalker,ndim,objective) #Machen Sie einen Sampler
sampler.run_mcmc(s0,5000) #Starten Sie die Probenahme

s_dist = sampler.flatchain #Holen Sie sich das Ergebnis der Probenahme
s = s_dist[sampler.flatlnprobability.argmax()]

y_pred_mean,y_pred_std = pred(x_train,y_train,x_pred,s)

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(x_true, func(x_true), label="True")
axes.scatter(x_train,y_train, label="Measured")
axes.plot(x_pred, y_pred_mean, label="Predict")
plt.fill_between(x_pred,y_pred_mean+y_pred_std,y_pred_mean-y_pred_std,facecolor="b",alpha=0.3)
axes.legend(loc="best")
axes.set_xlim(-2.0,12.0)

Es fühlt sich gut an.

3. Gaußsche Prozessregression unter Verwendung der GPy-Bibliothek

Es gibt eine praktische Bibliothek zur Gaußschen Prozessregression. Es wird empfohlen, da es eine Vielzahl von Kerneln enthält und leicht zu visualisieren ist.

import GPy
import GPy.kern as gp_kern

kern = gp_kern.RBF(input_dim=1)
gpy_model = GPy.models.GPRegression(X=x_train.reshape(-1, 1), Y=y_train.reshape(-1, 1), kernel=kern, normalizer=None)

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(x_true, func(x_true), c="k", label="True")
gpy_model.optimize()
gpy_model.plot(ax=axes)
axes.legend(loc="best")
axes.set_xlim(-2.0,12.0)

Zusammenfassung

Die Gaußsche Prozessregression ermöglichte die Schätzung einschließlich der Unsicherheit. In diesem Beitrag konnte ich nur den grundlegenden theoretischen Teil betrachten, aber die Gaußsche Prozessregression hat ein Problem mit dem Rechenvolumen, und es wurden verschiedene Diskussionen darüber geführt, wie es gelöst werden kann. Darüber hinaus gibt es viele Arten von Kernelfunktionen, nicht nur einen, und es ist erforderlich, sie entsprechend dem Analyseziel auszuwählen und entsprechend zu kombinieren.