Simulation des hinteren Teils eines Gaußschen Prozesses

8

Zum ersten Mal (entschuldigen Sie Ungenauigkeit / Fehler) habe ich mir Gaußsche Prozesse angesehen und mir dieses Video von Nando de Freitas genauer angesehen . Die Notizen sind hier online verfügbar .

Irgendwann zieht er 10 Zufallsstichproben aus einer multivariaten Normalen, die durch Konstruktion einer Kovarianzmatrix basierend auf einem Gaußschen Kern (Exponential der quadratischen Abstände in der x Achse) erzeugt wurden. Diese Zufallsstichproben bilden die vorherigen glatten Diagramme, die mit der Verfügbarkeit von Daten weniger verbreitet werden. Letztendlich besteht das Ziel darin, durch Modifizieren der Kovarianzmatrix und Erhalten der bedingten Gaußschen Verteilung an den interessierenden Punkten vorherzusagen.

Geben Sie hier die Bildbeschreibung ein

Der gesamte Code ist in einer ausgezeichneten Zusammenfassung von Katherine Bailey hier verfügbar , die wiederum ein Code-Repository von Nando de Freitas hier gutschreibt . Ich habe den Python-Code hier zur Vereinfachung veröffentlicht.

Es beginnt mit (anstelle von 10 oben) vorherigen Funktionen und führt einen "Abstimmungsparameter" ein.310

Ich habe den Code in Python und [R] übersetzt , einschließlich der Diagramme:

Hier ist der erste Codeabschnitt in [R] und die resultierende Darstellung von drei zufälligen Kurven, die über einen Gaußschen Kernel basierend auf der Nähe der Werte im Testsatz generiert wurden :x

! [Bildbeschreibung hier eingeben

Der zweite Teil des R-Codes ist haariger und beginnt mit der Simulation von vier Punkten von Trainingsdaten, die schließlich dazu beitragen, die Streuung zwischen den möglichen (vorherigen) Kurven um die Bereiche, in denen diese Trainingsdatenpunkte liegen, einzugrenzen. Die Simulation des Werts für diese Datenpunkte erfolgt als sin ( ) -Funktion. Wir können die "Verschärfung der Kurven um die Punkte" sehen:ySünde()

Geben Sie hier die Bildbeschreibung ein

Der dritte Teil des R-Codes befasst sich mit der Darstellung der Kurve der geschätzten Mittelwerte (das Äquivalent der Regressionskurve), die μ- Werten entspricht (siehe Berechnung unten), und ihrer Konfidenzintervalle:50 μ

Geben Sie hier die Bildbeschreibung ein

FRAGE: Ich möchte um eine Erklärung der Operationen bitten, die stattfinden, wenn ich vom vorherigen Hausarzt zum hinteren Hausarzt gehe.

Insbesondere möchte ich diesen Teil des R-Codes (im zweiten Block) verstehen, um die Mittel und SD zu erhalten:

# Apply the kernel function to our training points (5 points):

K_train = kernel(Xtrain, Xtrain, param)                          #[5 x 5] matrix

Ch_train = chol(K_train + 0.00005 * diag(length(Xtrain)))        #[5 x 5] matrix

# Compute the mean at our test points:

K_trte = kernel(Xtrain, Xtest, param)                            #[5 x 50] matrix
core = solve(Ch_train) %*% K_trte                                #[5 x 50] matrix
temp = solve(Ch_train) %*% ytrain                                #[5 x 1] matrix
mu = t(core) %*% temp                                            #[50 x 1] matrix

Es gibt zwei Kernel (einen von Zug ( ) gegen Zug ( a ), nennen wir es Σ a a mit seinem Cholesky ( ), L a a , der von nun an alle Cholesky orange färbt, und den zweiten von Zug ( a ) v - Test ( e ) , nennen sie es & Sgr; eine e ) und das geschätzte Mittel zur Erzeugung & mgr; für die 50 Punkte in dem Prüfgerät die Operation:eineinK_trainΣeineinCh_trainL.eineineineK_trteΣeineμ^50

(Gl. 1)μ^=[L.einein- -1[5×5]]Σeine[5×50]]]]T.L.einein- -1[5×5]]ytr[5×1]]Maße=[50×1]]
# Compute the standard deviation:

tempor = colSums(core^2)                                          #[50 x 1] matrix

# Notice that all.equal(diag(t(core) %*% core), colSums(core^2)) TRUE

s2 = diag(K_test) - tempor                                        #[50 x 1] matrix
stdv = sqrt(s2)                                                   #[50 x 1] matrix

(Gleichung 2)var^=diag(Σee)- -diag[[L.einein- -1[5×5]]Σeine[5×50]]]]T.[L.einein- -1[5×5]]Σeine[5×50]]]]]]=d[1111]]- -d[[L.einein- -1[5×5]]Σeine[5×50]]]]T.[L.einein- -1[5×5]]Σeine[5×50]]]]]]Maße=[50×1]]

Wie funktioniert das?

μ^

Ch_post_gener = chol(K_test + 1e-6 * diag(n) - (t(core) %*% core))
m_prime = matrix(rnorm(n * 3), ncol = 3)
sam = Ch_post_gener %*% m_prime
f_post = as.vector(mu) + sam 

(Gleichung 3)fPost=μ^+[L.ee[50×50]]- -[[L.einein- -1[5×5]]Σeine[5×50]]]]T.[L.einein- -1[5×5]]Σeine[5×50]]]]]]]][N.(0,1)[50×3]]]]Maße=[50×3]]
Antoni Parellada
quelle
Sollten die Konfidenzintervalle in der letzten Darstellung nicht an den bekannten Punkten "kneifen"?
GeoMatt22
@ GeoMatt22 Sie tun es irgendwie, meinst du nicht auch?
Antoni Parellada

Antworten:

2

eeineine

[eine]]N.([μeinμe]],[ΣeineinΣeineΣeineT.Σee]])

E.(x1|x2)=μ1+Σ12Σ22- -1(x2- -μ2)[50×50]]Σeinein[50×5]]Σeinewird eine transponierte notwendig sein, um die Matrizen kongruent zu machen in:

E.(e|ein)=μe+ΣeineT.Σeinein- -1(y- -μein)
μein=μe=0

E.(e|ein)=ΣeineT.Σeinein- -1ytr

Geben Sie die Cholesky-Zerlegung ein (die ich wieder wie in OP in Orange codieren werde):

E.(e|ein)=ΣeineT.Σeinein- -1ytr<- -- -α- -- ->=ΣeineT.(L.eineinL.eineinT.)- -1ytr=ΣeineT.L.einein- -T.L.einein- -1ytr(*)=ΣeineT.L.einein- -T.L.einein- -1ytr<- -m- ->

m=L.einein- -1ytrL.eineinm=ytrm


Geben Sie hier die Bildbeschreibung ein


B.T.EINT.=(EINB.)T.

μ^=[L.einein- -1[5×5]]Σeine[5×50]]]]T.L.einein- -1[5×5]]ytr[5×1]]=(ΣeineT.L.einein- -T.)(L.einein- -1ytr)Maße=[50×1]]

angesichts dessen

[L.einein- -1[5×5]]Σeine[5×50]]]]T.=Σeine[50×5]]T.L.einein- -1T.[5×5]]

Eine ähnliche Argumentation würde auf die Varianz angewendet, beginnend mit der Formel für die bedingte Varianz in einem multivariaten Gaußschen:

veinr(x1|x2)=Σ11- -Σ12Σ22- -1Σ21

was in unserem Fall wäre:

varμ^e=Σee- -ΣeineT.Σeinein- -1Σeine=Σee- -ΣeineT.[L.eineinL.eineinT.]]- -1Σeine=Σee- -ΣeineT.[L.einein- -1]]T.L.einein- -1Σeine=Σee- -[L.einein- -1Σeine]]T.L.einein- -1Σeine

und Erreichen von Gleichung (2):

varμ^e=d[K.ee- -[L.einein- -1[5×5]]Σeine[5×50]]]]T.[L.einein- -1[5×5]]Σeine[5×50]]]]]]Maße=[50×1]]

Wir können sehen, dass Gleichung (3) im OP eine Möglichkeit ist, hintere zufällige Kurven zu erzeugen, die von den Daten abhängig sind (Trainingssatz), und eine Cholesky-Form zu verwenden, um drei multivariate normale zufällige Ziehungen zu erzeugen :

fPost=μ^+[varμ^e]][rnorm(0,1)]]=μ^+[L.ee[50×50]]- -[[L.einein- -1[5×5]]Σeine[5×50]]]]T.[L.einein- -1[5×5]]Σeine[5×50]]]]]]]][rand.norm's[50×3]]]]Maße=[50×3]]
Antoni Parellada
quelle
1
Ist das aus einem Buch oder Papier? Haben Sie eine robuste Methode, um den bedingten Mittelwert und die Varianz zu berechnen, wenn die Kovarianzmatrix EXTREM schlecht konditioniert ist (ohne jedoch nahezu abhängige (nahe gelegene) Datenpunkte zu löschen oder zusammenzuführen), und zwar mit doppelter Genauigkeit? Die Mehrfachpräzision in Software funktioniert, hat jedoch eine Verlangsamung um 2,5 bis 3 Größenordnungen im Vergleich zur Hardware-Doppelpräzision, sodass selbst ein "langsamer" Algorithmus mit doppelter Genauigkeit gut ist. Ich glaube nicht, dass Cholesky es schneidet. Ich denke auch nicht, dass QR dies tut, wenn die Kovarianzmatrix sehr schlecht konditioniert ist. Bei Verwendung von Standard-Backsolves scheint eine mehrfache Präzision erforderlich zu sein.
Mark L. Stone