Gaußsches Prozessvorhersageintervall

Ich werde Ihre Frage im Rahmen von Bayes beantworten. Wenn Sie speziell eine frequentistische Lösung benötigen, können Sie eine erhalten, indem Sie meine Antwort leicht ändern. Ich denke jedoch, dass dies die tatsächliche Unsicherheit unterschätzen wird: Sie würden einen vollständig frequentistischen Ansatz benötigen, aber ich weiß nicht, wie ich das in diesem speziellen Fall tun soll .

Um das Bayesian GPR-Framework (Gaussian Process Regression) kurz zusammenzufassen, nehmen Sie das Modell an

y = f (x | θ) + ϵ

$y=f(x|\boldsymbol{\theta})+\epsilon$

wobei , dh Die latenten Variablen oder Funktionswerte werden als Gaußscher Prozess unter der Bedingung der Hyperparameter , und ist das übliche iid-Gaußsche Rauschen. $f(x|\boldsymbol{\theta})\sim \mathcal{GP}(\mu(x|\boldsymbol{\theta}),k(x,x'|\boldsymbol{\theta}))$ $\boldsymbol{\theta}$ $\epsilon\sim\mathcal{N}(0,\sigma^2)$

Eigentlich ist ein Hyperparameter, also gehört er wirklich zu , aber ich wollte unterstreichen, dass GPR normalerweise eine triviale Kovarianzstruktur für das Rauschen annimmt. $\sigma^2$ $\boldsymbol{\theta}$

Die posteriore prädiktive Verteilung von an einem neuen Punkt , bedingt durch Daten und auf Hyperparametern ist . Nehmen wir nun an, dass die mittlere Funktion des Gaußschen Prozesses Null ist: Der allgemeine Fall kann ebenfalls behandelt werden, aber versuchen wir, die Dinge einfach zu halten. Dann bekommen wir mit der üblichen GPR-Maschinerie $y_*$ $x_*$ $\{(x_1,y_1,)\dots,(x_d,y_d)\}=(\mathbf{x},\mathbf{y})$ $\boldsymbol{\theta}$ $p(y_*|\boldsymbol{\theta},\mathbf{y})$

p (f_{*} | θ, y) = N (k_{*}^{T} (K + σ^{2} I)^{- 1} y, k (x_{*}, x_{*}) - k_{*}^{T} (K + σ^{2} I)^{- 1} k_{*})

$p(f_*|\boldsymbol{\theta},\mathbf{y}) = \mathcal{N}(\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{y},k(x_*,x_*)-\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{k}_*)$

K = (\begin{matrix} k (x_{1}, x_{1}; θ) & k (x_{1}, x_{d}; θ) \\ ⋱ \\ k (x_{1}, x_{d}; θ) & k (x_{d}, x_{d}; θ) \end{matrix})

$K=\pmatrix{k(x_{1},x_{1};\boldsymbol{\theta})& & k(x_{1},x_{d};\boldsymbol{\theta}) \\ & \ddots & \\k(x_{1},x_{d};\boldsymbol{\theta})& & k(x_{d},x_{d};\boldsymbol{\theta}) }$

k_{*} = (\begin{matrix} k (x_{*}, x_{1}; θ) \\ ⋮ \\ k (x_{*}, x_{d}; θ) \end{matrix})

$\mathbf{k}_*=\pmatrix{k(x_*,x_{1};\boldsymbol{\theta}) \\ \vdots \\k(x_*,x_{d};\boldsymbol{\theta})}$

dh abhängig von beobachteten Daten und Hyperparametern ist die Verteilung der latenten Variablen an einem neuen Punkt immer noch Gaußsch, wobei der Mittelwert und die Standardabweichung oben gezeigt sind.

Wir sind jedoch an der Verteilung einer neuen Beobachtung interessiert , nicht einer neuen latenten Variablen. Dies ist einfach, da in unserem Modell das Rauschen additiv ist, unabhängig von allen anderen Variablen und normalerweise mit dem Mittelwert Null und der Varianz verteilt ist. Daher müssen wir nur die Rauschvarianz addieren: $y_*$ $\sigma^2$

p (y_{*} | θ, y) = N (k_{*}^{T} (K + σ^{2} I)^{- 1} y, k (x_{*}, x_{*}) - k_{*}^{T} (K + σ^{2} I)^{- 1} k_{*} + σ^{2})

$p(y_*|\boldsymbol{\theta},\mathbf{y}) = \mathcal{N}(\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{y},k(x_*,x_*)-\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{k}_*+\sigma^2)$

Beachten Sie, dass ich eine einzelne neue Beobachtung Betracht , sodass die Verteilung nur ein univariater Gaußscher Wert ist und die Varianz tatsächlich eine Varianz und keine a ist Varianz-Kovarianz-Matrix. $y_*$ $p(y_*|\boldsymbol{\theta},\mathbf{y})$

Um diesen Ausdruck tatsächlich verwenden zu können, benötigen Sie Werte für die Hyperparameter, die nicht bekannt sind. Hieraus gibt es zwei Möglichkeiten:

(die häufigste Lösung) Die Hyperparameter werden durch MLE oder MAP geschätzt, und der obige Ausdruck wird verwendet. Dieser Ansatz vernachlässigt die Unsicherheit bei der Schätzung der Hyperparameter vollständig und scheint daher nicht sehr sicher zu sein.
In einem vollständig Bayes'schen Ansatz interessieren Sie sich nicht wirklich für , sondern für die prädiktive Verteilung von bei , die durch erhalten wird nach Integration der Hyperparameter: $p(y_*|\boldsymbol{\theta},\mathbf{y})$ $y_*$ $\mathbf{y}$ $p(y_*|\boldsymbol{\theta},\mathbf{y})$

$p (y_{*} | y) = \int p (y_{*}, θ | y) d θ = \int p (y_{*} | θ, y) p (θ | y) d θ$ $p(y_*|\mathbf{y})=\int{p(y_*,\boldsymbol{\theta}|\mathbf{y})}d\boldsymbol{\theta}=\int{p(y_*|\boldsymbol{\theta},\mathbf{y})p(\boldsymbol{\theta}|\mathbf{y})}d\boldsymbol{\theta}$

Hier gibt es zwei Probleme: Bei einer vorherigen Verteilung für die Hyperparameter ist die posteriore Verteilung , die im Integral erscheint, gegeben nicht bekannt, muss aber unter Verwendung des Bayes-Theorems abgeleitet werden, was für die meisten Hyperprioren bedeutet, dass ein MCMC ausgeführt werden muss. Daher haben wir keinen expliziten Ausdruck für , sondern nur Beispiele aus der MCMC. Und selbst wenn wir einen Ausdruck für , wäre es unmöglich , das Integral in geschlossener Form in zu bewerten meiste Fälle. Die Lösung ist eine hierarchische Bayes-Simulation: für jede Probe $p(\boldsymbol{\theta})$ $p(\boldsymbol{\theta}|\mathbf{y})$ $p(\boldsymbol{\theta}|\mathbf{y})$ $p(\boldsymbol{\theta}|\mathbf{y})$ $p(y_*|\mathbf{y})$ $\hat{\boldsymbol{\theta}}_i$ aus mit dem MCMC erhalten Sie eine Stichprobe aus . Verwenden Sie diese Stichproben , um ein HPD-Intervall für zu schätzen , und schon sind Sie da. $p(\boldsymbol{\theta}|\mathbf{y})$ $y^*_i$ $p(y_*|\hat{\boldsymbol{\theta}}_i,\mathbf{y})$ $m$ $y^*_i$ $y_*$

Aus intuitiver Sicht werden bei der zweiten Lösung Stichproben aus einer Verteilung gezogen, bei der die Hyperparameter "nicht festgelegt" sind, aber entsprechend ihrer posterioren Verteilung zufällig variieren dürfen . Somit berücksichtigt das im zweiten Fall erhaltene Vorhersageintervall die Unsicherheit aufgrund unseres Mangels an Wissen über die Hyperparameter. $p(\boldsymbol{\theta}|\mathbf{y})$

DeltaIV
quelle

Tut mir leid, aber bitte können Sie in der Antwort kurz erläutern, wie die Vorhersageformel erhalten wird, wie Sie es in den Kommentaren zur anderen Antwort getan haben. Dadurch wird die Antwort in sich geschlossen und für jeden nützlich, der sie in Zukunft sucht. Tolle Diskussion über Hyperparameter-Probleme übrigens! Danke :)

EHH

@EHH kein Problem, ich habe das Bit hinzugefügt, auf das Sie sich bezogen haben.

DeltaIV

@DeltalV Super, das hat mir wirklich geholfen, einige Dinge zu klären, über die ich mich gewundert habe! Vielen Dank!

EHH

@ Mathews24 das ist eine andere Frage und die Lebenslaufrichtlinie ist eine Frage pro Beitrag. Durchsuchen Sie die Website, um festzustellen, ob eine solche Frage bereits gestellt wurde. Andernfalls können Sie selbst eine neue Frage stellen.

DeltaIV

@ DeltaIV Sie geben an, "wir sind an der Verteilung einer neuen Beobachtung interessiert ". Ich nehme an, es ist kontextabhängig, aber gibt es allgemeine Regeln, wenn man sich für versus interessiert ? Was ist die physikalische Interpretation für ? Wenn beispielsweise Messungen aus einer Diagnose mit einem Fehler , würde uns die Modellierung von tatsächlich das Modell für das zugrunde liegende physikalische Phänomen geben? Warum würden wir in diesem Fall bevorzugen ?

y_{*}

$y_∗$

y_{*}

$y_*$

f_{*}

$f_*$

f_{*}

$f_*$

y

$y$

σ

$\sigma$

f_{*}

$f_*$

y_{*}

$y_*$

Mathews24

Wenn Sie sich auf die Bayes'sche Regression mit Gaußscher Wahrscheinlichkeit beziehen, ist die hintere Verteilung eines Gaußschen Prozesses Gauß'sch: wobei die und die Datenwerte sind und und mit Bayes'scher Inferenz berechnet werden: wobei der Kernelvektor zwischen und und mit

p (f (x) ∣ X_{n}, Y_{n}) = N (μ_{n} (x), σ_{n}^{2} (x)),

$p(f(x)\mid X_n,Y_n) = \mathcal{N}\big(\mu_n(x),\sigma_n^2(x)\big)\,,$

X_{n}

$X_n$

Y_{n}

$Y_n$

μ_{n}

$\mu_n$

σ_{n}^{2}

$\sigma_n^2$

μ_{n} (x) = k_{n} (x)^{⊤} C_{n}^{- 1} Y_{n} and σ_{n}^{2} (x) = k (x, x) - k_{n} (x)^{⊤} C_{n}^{- 1} k_{n} (x),

$\mu_n(x) = \mathbf{k}_n(x)^\top C_n^{-1}Y_n \text{ and } \sigma_n^2(x) = k(x,x) - \mathbf{k}_n(x)^\top C_n^{-1} \mathbf{k}_n(x)\,,$

k_{n} (x) = [k (x_{t}, x)]_{x_{t} \in X_{n}}

$\mathbf{k}_n(x) = [k(x_t, x)]_{x_t \in X_n}$

x

$x$

X_{n}

$X_n$

C_{n} = K_{n} + η^{2} I

$C_n = K_n + \eta^2 I$

η^{2}

$\eta^2$ die Standardabweichung des Beobachtungsrauschens und der Kernelmatrix (siehe das zweite Kapitel von Rasmussen und Williams 'Buch ).

K_{n} = [k (x_{t}, x_{t^{'}})]_{x_{t}, x_{t^{'}} \in X_{n}}

$K_n=[k(x_t,x_{t'})]_{x_t,x_{t'} \in X_n}$

Daher ist ein ~ 95% -Konfidenzintervall für einfach . $x$ $\mu_n(x) \pm 2 \sigma_n(x)$

Emile
quelle

Es hört sich so an, als würden Sie annehmen, dass (a) das OP ein Bayes'sches Vorhersageintervall wünscht und (b) bereit ist, vor dem Prozess ein normales Konjugat zu übernehmen.

whuber

Tatsächlich. Ich nehme diese Annahme in meine Antwort auf.

Emile

Wie kann U σ2n (x) finden? Kennst du die Gleichung dafür? Vielen Dank

Wis

Diese Antwort ist absolut falsch. Das OP sagte, dass sie wissen, wie man die Konfidenzintervalle erhält, aber die Vorhersageintervalle erhalten möchten, bei denen es sich um die 95% -Wahrscheinlichkeitsintervalle für zukünftige Beobachtungen handelt . Ich empfehle @ raK1, dies als akzeptierte Antwort zu entfernen.

w i t h n o i s e

$\mathbf{with ~noise}$

EHH

@DeltaIV Ich glaube, dass die zukünftige Messvarianz, die die Summe aus Funktion und Rauschvarianz ist, die Antwort ist, nach der das OP gesucht hat (ich dachte, es könnte dies sein, hatte aber keine Referenz dafür gefunden und hoffte daher auf jemanden, der sich sicher war zu beantworten), also sollten Sie dies vielleicht als Antwort auf die Frage schreiben. In Bezug auf Ihren anderen Punkt stimme ich voll und ganz zu, dass die Annahme, dass Hyperparameter bekannt sind, die Unsicherheit unterschätzt. Pragmatisch gesehen wird dies jedoch bei der praktischen Anwendung von Hausärzten getan, sodass ich in diesem Zusammenhang geantwortet habe.

EHH

Gaußsches Prozessvorhersageintervall

Antworten: