[PYTHON] Ich habe GP mit Numpy geschrieben

Einführung

Schön, Sie kennenzulernen. Mein Name ist Mimura und ich werde für den 21. Tag des Adventskalenders der NTT Docomo Service Innovation Department verantwortlich sein. Ich werde dieses Jahr als Freiwilliger von Docomo einen Adventskalender schreiben, daher möchte ich einen Artikel schreiben.

Kennen Sie plötzlich das Regressionsproblem? Das Regressionsproblem wird an verschiedenen Stellen verwendet, beispielsweise bei der Prognose von Aktienkursen. DoCoMo befasst sich auch an verschiedenen Stellen mit diesem Regressionsproblem.

In diesem Artikel ist die Gaußsche Prozessregression eine der Methoden zur Lösung dieses Regressionsproblems, und ich werde sie implementieren. Der Zweck dieser Zeit ist es, sie mit Numpy zu implementieren, nicht mit einer theoretischen Erklärung. Danke. Für diejenigen, die die spezifischen Inhalte dieses Bereichs kennen möchten ["Gauß-Prozess und maschinelles Lernen (Professional Series für maschinelles Lernen)"](https://www.amazon.co.jp/%E3%82%AC%E3%82% A6% E3% 82% B9% E9% 81% 8E% E7% A8% 8B% E3% 81% A8% E6% A9% 9F% E6% A2% B0% E5% AD% A6% E7% BF% 92- % E6% A9% 9F% E6% A2% B0% E5% AD% A6% E7% BF% 92% E3% 83% 97% E3% 83% AD% E3% 83% 95% E3% 82% A7% E3 % 83% 83% E3% 82% B7% E3% 83% A7% E3% 83% 8A% E3% 83% AB% E3% 82% B7% E3% 83% AA% E3% 83% BC% E3% 82 % BA-% E6% 8C% 81% E6% A9% 8B-% E5% A4% A7% E5% 9C% B0 / dp / 4061529269 / ref = sr_1_1? __Mk_ja_JP =% E3% 82% AB% E3% 82% BF% E3% 82% AB% E3% 83% 8A & Schlüsselwörter =% E3% 82% AC% E3% 82% A6% E3% 82% B9% E9% 81% 8E% E7% A8% 8B & qid = 1575191954 & sr = 8-1 ](Https://www.amazon.co.jp/ Gauß-Prozess und maschinelles Lernen - Maschinelles Lernen Professionelle Serie-Mochihashi-Erde / dp / 4061529269 / ref = sr_1_1? __Mk_ja_JP = Katakana & keywords = Gauß-Prozess & qid = 1575191954 & sr = 8 -1 "Aufruf" Gaußscher Prozess und maschinelles Lernen (Machine Learning Professional Series) "")

Wenn Sie tatsächlich den Gaußschen Prozess verwenden, verwenden Sie auch GPyTorch oder GPy! Schließlich konnten wir diesmal weder die Kernelfunktionen optimieren noch die Parameter anpassen. Es tut uns leid…

Eine kurze Erklärung zur Lösung des Regressionsproblems

Einfach ausgedrückt lautet das einfache Regressionsproblem "Ich möchte einige Daten $ Y $ aus einigen beobachtbaren Daten $ X $ vorhersagen". Eine Möglichkeit, dieses Problem zu lösen, besteht darin, eine Funktion wie $ y = Ax + B $ vorzubereiten. Wenn wir die Variable $ x_i $ beobachten können, können wir $ y_i $ vorhersagen.

In der Realität ist es jedoch schwierig, es mit einer einfachen Funktion wie $ y = Ax + B $ auszudrücken. Um dies zu lösen, machen Sie die Formel schwierig und verbessern Sie die Ausdruckskraft wie $ y = Ax ^ 2 + Bx + C $.

Dieses Mal werden wir diese komplizierte Funktion als $ y = W ^ T \ phi (x) $ ausdrücken. Wenn zum Beispiel $ W = [A, B, C], \ phi (x) = [x ^ 2, x ^ 1, x ^ 0] $, dann wird $ y = Ax ^ 2 + Bx + C $ gut ausgedrückt. tun können. Mit anderen Worten, wenn Sie $ \ phi (x) $ schwierig machen und $ W $ dafür finden können, können Sie $ Y $ gut vorhersagen! !! !! !!

Yay! !! !! !! !! !! !! !! !! !! !! !! !! !! !! !! !! !! !!

…

Wenn es hier endet, wird es gegen das Thema sein, also werde ich den Gauß-Prozess von hier aus erklären.

Eine kurze Beschreibung des Gaußschen Prozesses

Betrachten Sie zunächst, $ \ phi (x) $ von $ y = W ^ T \ phi (x) $ zu einer Gaußschen Verteilung zu machen. Mit anderen Worten, versuchen Sie $ \ phi (x) = \ exp (- \ frac {(x- \ mu_h) ^ 2} {\ sigma}) $. Das Regressionsproblem kann noch gelöst werden.

Bei dieser Methode wird der Rechenaufwand jedoch aufgrund des Fluches der Dimensionen enorm. Um dieses Problem zu lösen, sollten Sie eine Methode in Betracht ziehen, um den erwarteten Wert in $ W $ zu übernehmen und ihn zu integrieren und aus dem Modell zu entfernen.

Wenn $ w $ aus einer Gaußschen Verteilung mit einem Durchschnitt von $ 0 $ und einer Varianz von $ \ lambda ^ 2 I $ als $ y = \ Phi W $ erzeugt wird, ist es $ w \ sim N (0, \ lambda ^ 2 I) $. .. Dies bedeutet, dass $ y $ "eine lineare Transformation von $ W $ ist, die einer Gaußschen Verteilung mit der konstanten Matrix $ \ Phi $ folgt".

Zu diesem Zeitpunkt sind der erwartete Wert und die Kovarianz

• E[y]=E[\Phi W]=\Phi E[ww^T]=0

• \Sigma = E[yy^T]-E[y]E[y]^T=E[(\Phi w)(\Phi w)^T]=\Phi E[ww^T]\Phi^T=\lambda^2\Phi\Phi^T

Wird sein

Aus diesem Ergebnis folgt $ y $ einer multivariaten Gaußschen Verteilung von $ y \ sim N (0, \ lambda ^ 2 \ Phi \ Phi ^ 2) $.

$ K = \ lambda ^ 2 \ Phi \ Phi ^ T $ und $ K_ {n, n '} = \ lambda ^ 2 \ phi (X_n) \ phi (X_ {n'}) $ und der Durchschnitt von $ y $ Wenn Sie auf $ 0 $ normalisieren, erhalten Sie $ y \ sim N (0, K) $. Trotzdem ist die Berechnung von $ K_ {n, n '} = \ phi (X_n) \ phi (X_ {n'}) $ immer noch schwer ... Aber hier ist das berühmte "Lass uns Kernel-Tricks anwenden!" Die Idee ist, dass $ K_ {n, n '} $ berechnet werden kann, ohne hart arbeiten zu müssen, um $ \ phi (X_n) $ zu berechnen.

Wenn Sie zu diesem Punkt kommen, können Sie mehr Daten erhalten

D=\{(x_1,y_1),(x_2,y_2),\cdots,(x_N,y_N)\} \\

Was passiert mit $ y ^ {new} $, wenn neue Daten beobachtet werden? $ X ^ {new} $ = y ^ {new} = w ^ T \ phi (x ^ {new) }) Alles was Sie tun müssen, ist $ zu lösen.

Um dies zu erreichen, ist $ y '= (y_1, y_2, \ cdots, y_N, y ^ {new}), X = (x_1, x_2, \ cdots, x_N, x ^ {new}) $ Es wird y \ sim N (0, K ') $.

Diesen Weg

\begin{pmatrix}
y \\ 
y^{new}
\end{pmatrix}
\sim 
N\left(0,\begin{pmatrix}
K,k_{new} \\ 
k_{new}^T,k_{new,new}
\end{pmatrix}\right)

Kann geschrieben werden als. Hier sind $ k_ {new} und k_ {new, new} $ wie folgt.

\left\{ 
\begin{array}{}
k_{new} &=(k(x^{new},x_1),k(x^{new},x_2),\cdots,k(x^{new},x_N))^T \\
k_{new,new}&= (k(x^{new},x^{new}))
\end{array}
\right.

Hier

\begin{pmatrix}
y_1 \\ 
y_2
\end{pmatrix}
\sim N\left(
\begin{pmatrix}
\mu_1 \\ 
\mu_2
\end{pmatrix}
,
\begin{pmatrix}
\Sigma_{1,1},\Sigma_{1,2} \\
\Sigma_{2,1},\Sigma_{2,2}
\end{pmatrix}
\right)

Als es die Formel gabp(y_2|y_1)Zup(y_2|y_1)=N(\mu_2+\Sigma_{2,1}\Sigma_{1,1}^{-1}(y_1-\mu_1),\Sigma_{2,2}-\Sigma_{2,1}\Sigma_{2,1}^{-1}\Sigma_{1,2})Es ist bekannt, dass es ausgedrückt werden kann durch. Diesmal\mu_1=0,\mu_2=0Damitp(y_2|y_1)=N(\Sigma_{2,1}\Sigma_{1,1}^{-1}y_1,\Sigma_{2,2}-\Sigma_{2,1}\Sigma_{2,1}^{-1}\Sigma_{1,2})Es kann ausgedrückt werden durch.

Es war so lange, aber mit anderen Worten

p(y^{new}|X^{new},D)=N(k_{new}^TK^{-1}y,k_{new,new}-k_{new}^TK^{-1}k_{new})

Alles was Sie tun müssen, ist zu implementieren! !! !!

Lassen Sie es uns tatsächlich implementieren!

p(y^{new}|X^{new},D)=N(k_{new}^TK^{-1}y,k_{new,new}-k_{new}^TK^{-1}k_{new})

Ich habe gerade herausgefunden, dass ich dies implementieren sollte.

Auch $ K, k_ {new}, k_ {new, new} $

\left\{ 
\begin{array}{}
K   &=\left[\begin{array}{}
k(x_1,x_1),k(x_1,x_2),\cdots,k(x_1,x_N) \\
k(x_2,x_1),k(x_2,x_2),\cdots,k(x_2,x_N) \\
\cdots \\
k(x_N,x_1),k(x_N,x_2),\cdots,k(x_N,x_N)
\end{array}
\right] \\
k_{new} &=(k(x^{new},x_1),k(x^{new},x_2),\cdots,k(x^{new},x_N))^T \\
k_{new,new}&= (k(x^{new},x^{new}))
\end{array}
\right.

war.

k(x,x')=\theta_1\exp \left(-\frac{(x-x')^2}{\theta_2}\right)+\theta_3\delta(x,x')

Dann ist die Kernelfunktion

def RBF(X,X_):
    theta_1  = 1
    theta_2  = 0.2
    theta_3  = 0.01
    distance =  ((X - X_)**2).sum(-1)
    if distance.shape[0] == distance.shape[1]:
        return theta_1 * np.exp(-1 * distance/theta_2) + theta_3 * np.eye(len(X))[:,:X.shape[1]]
    else:
        return theta_1 * np.exp(-1 * distance/theta_2) 

Sie können mit schreiben.

Verwenden wir dies, um $ K ^ {-1}, k_ {new}, k_ {new, new} $ zu berechnen. Berechnen Sie zunächst $ K ^ {-1} $.

X_    = np.array([X for _ in range(len(X))])
K     = RBF(X_,X_.transpose(1,0,2)) 
inv_K = np.linalg.inv(K)

Berechnen Sie dann $ k_ {new} $ und $ k_ {new, new} $.

X__        = np.array([X      for _ in range(len(Test_X))])
Test_X__   = np.array([Test_X for _ in range(len(X))]).transpose(1,0,2)
Test_X__a  = np.array([Test_X for _ in range(len(Test_X))]).transpose(1,0,2)
k          = RBF(X__,Test_X__)
k__        = RBF(Test_X__a,Test_X__a)

Schließlich

p(y^{new}|X^{new},D)=N(k_{new}^TK^{-1}y,k_{new,new}-k_{new}^TK^{-1}k_{new})

Zu generieren aus

Y_result   = np.random.multivariate_normal(k.dot(inv_K).dot(Y),k__ - k.dot(inv_K).dot(k.T))

Ruf einfach an

Experimentieren Sie in New York mit Daten zur Fahrradfreigabe

Verwenden Sie die New Yorker Bike Share-Daten (https://www.citibikenyc.com/system-data), um Ihr Bild vorherzusagen. Dieses Mal habe ich versucht, die Anzahl der zurückgegebenen Fahrräder am 30. Juni von 8.00 bis 12.00 Uhr zu visualisieren.

Aus Gründen der Übersichtlichkeit wird es logarithmisch auf die zurückgegebene Zahl normiert.

Lassen Sie es uns mit dem Gaußschen Prozess visualisieren Dann sieht es so aus.

Der Stern ist die Position des Hafens. $ 0 $ und darunter werden für eine einfache Anzeige mit $ 0 $ ergänzt.

Die Ergebnisse zeigen, dass im Zentrum eine hohe Nachfrage besteht. Es wurde gefolgert, dass diesmal eine hohe Nachfrage bestehen würde, da aufgrund der großen Streuung verschiedene Werte an Positionen ohne Ports abgetastet werden. Das Gute am Gaußschen Prozess ist, dass Sie nicht nur vorhersagen können, wo keine Daten vorhanden sind, sondern auch die Varianz dieses Teils berechnen können.

Schließlich

Sie können das Regressionsproblem wie folgt lösen. Warum versuchen Sie es nicht einmal?