Was ist Bayes'sche Lineare Regression? (Teil 1)

Matthias Werner

Bayes' Theorem

Bayes'sche Regressionsmethoden sind sehr leistungsfähig, da sie uns nicht nur Punktschätzungen von Regressionsparametern geben, sondern vielmehr eine vollständige Verteilung über diese Parameter liefern. Das kann so verstanden werden, dass man nicht nur ein Modell, sondern eine ganze Familie von Modellen lernt und ihnen je nach ihrer Wahrscheinlichkeit, korrekt zu sein, unterschiedliche Gewichte zuweist. Da diese Gewichtsverteilung von den beobachteten Daten abhängt, können Bayes'sche Methoden uns eine Unsicherheitsquantifizierung unserer Vorhersagen geben, die darstellt, was das Modell aus den Daten lernen konnte. Das Unsicherheitsmaß könnte z.B. die Standardabweichung der Vorhersagen aller Modelle sein. Dies ist etwas, was Punktschätzer nicht standardmäßig liefern. Zu wissen, was das Modell nicht weiß, trägt dazu bei, die KI besser erklärbar zu machen. Um die Grundidee der Bayes'schen Regression zu klären, bleiben wir bei der Diskussion der Bayes'schen linearen Regression (BLR).

BLR ist der Bayes'sche Ansatz zur linearen Regressionsanalyse. Wir werden mit einem Beispiel beginnen, um die Methode zu motivieren. Zur Verdeutlichung stellen wir dann einige nicht-bayesianische Methoden vor, mit denen der Leser vielleicht bereits vertraut ist, und diskutieren, wie sie sich zur Bayes'schen Regression verhalten.

Im Folgenden gehe ich davon aus, dass Sie über elementare Kenntnisse der linearen Algebra und der Stochastik verfügen.

Fangen wir also an!

Ein einführendes Beispiel

Die Aufgabe, die wir lösen wollen, ist die folgende: Angenommen, wir erhalten einen Satz $$D$$ von $$N_D$$ unabhängig generierten, rauschbehafteten Datenpunkten $$D := \{ (x_i, y_i)\}{i=1}^{N_D}$$, wobei $$x_i \in \mathbb{R}^m$$ und $$y_i \in \mathbb{R}^n$$ (in unserem Beispiel $$m=n=1$$). Wir wissen, dass es eine gewisse Beziehung zwischen den $$x_i$$ und den $$y_i$$ gibt, aber wir wissen nicht, was diese Beziehung ist. Wir werden eine vereinfachende Annahme treffen, nämlich dass die Beziehung linear ist, aber wir haben auch ein gewisses Messrauschen vorliegen. Formal können wir dies so schreiben:

$$\begin{equation} \begin{aligned} y_i &= w_1 x_i + w_0 + \epsilon_i = \vec{w}^T \vec{x}_i + \epsilon_i \\ \epsilon_i &\sim \mathcal{N}(0, \sigma_\epsilon^2 ),  \end{aligned} \label{eqModel} \end{equation}$$

wobei wir der Einfachheit halber $$\vec{x}_i := (1, x_i)^T$$ und $$\vec{w} := (w_0, w_1)^T$$ definieren. $$w_0$$ und $$w_1$$ sind die Modellparameter, die wir schätzen wollen. Alle $$\epsilon_i$$ werden unabhängig voneinander aus einer Normalverteilung mit der Standardabweichung $$\sigma_\epsilon$$ gezogen. Wir wollen das Rauschen $$\epsilon_i$$ nicht modellieren, sondern es als Ausrede dafür benutzen, dass unser Modell die Daten nicht perfekt repräsentiert.

Nehmen wir an, wir kennen bereits die posteriore Verteilung $$p(\vec{w}|D)$$, die das kodiert, was unserer Meinung nach die Parameter $$\vec{w}$$ sein könnten, nachdem wir die Daten $$D$$ beobachtet haben (wir werden in einer Minute erfahren, wie wir sie erhalten können). Für eine gegebene Eingabe $$\vec{x}$$ erhalten wir eine gesamte Verteilung über die Vorhersagen $$y = \vec{w}^T \vec{x}$$ gewichtet mit $$p(\vec{w}|D)$$. Diese Verteilung wird als prädiktive Verteilung bezeichnet.

Wir können die prädiktive Verteilung visualisieren, indem wir den Mittelwert und die Standardabweichung der Vorhersage für einen gegebenen $$x$$ aufzeichnen. Das Schöne an der BLR ist, dass wir diese Momente analytisch berechnen können. Für die mittlere Vorhersage erhalten wir

$$\begin{equation} \begin{aligned} \mu_{y|\vec{x}} &= \int d\vec{w} \ p(\vec{w}|D) \ \vec{w}^T \vec{x} \\ &= \vec{\mu}_w \vec{x}, \end{aligned} \end{equation}$$

wo wir die Linearität sowohl des Erwartungswertes als auch des Skalarproduktes verwenden. Für die Varianz müssen wir Folgendes berechnen:

$$\begin{equation} \begin{aligned} \sigma_{y|\vec{x}}^2 &= \int d\vec{w} \ p(\vec{w}|D) \ \left( \vec{w}^T \vec{x} - \vec{\mu}_w^T \vec{x} \right)^2 \\ &= \vec{x}^T \left[ \int d\vec{w} \ p(\vec{w}|D) \ \left( \vec{w} \vec{w}^T - \vec{\mu}_w \vec{\mu}_w^T \right) \right] \vec{x} \\ &= \vec{x}^T \Sigma_w \vec{x}, \end{aligned} \end{equation}$$

wobei $$\vec{\mu}_w$$ und $$\Sigma_w$$ der mittlere Vektor und die Kovarianzmatrix von $$p(\vec{w}|D)$$ sind.

Darüber hinaus können wir auch eine Stichprobe von $$p(\vec{w}|D)$$ nehmen und die resultierenden Modelle darstellen. Da $$\Sigma_w$$ positiv-definit ist (wir nehmen an, dass dies wahr ist, ohne es hier zu beweisen, aber es hat etwas damit zu tun, dass es die Summe einer positiven semi-definiten Matrix und einer identischen Matrix ist), können wir es so zerlegen, dass

$$\begin{equation} \Sigma_w = GG^T. \end{equation}$$

Wenn wir einen Vektor $$\vec{z}$$ aus einer multivariaten Standardnormalverteilung ziehen, d.h.

$$\begin{equation} \begin{aligned} \vec{z} &= [z_1, z_2, ... ]^T \\ z_i &\sim \mathcal{N}(0, 1) \end{aligned} \end{equation}$$

können wir es in eine Probe des Posteriors umwandeln:

$$\begin{equation} \vec{w} = G^T \vec{z} + \vec{\mu}_w \end{equation}.$$

Before observing data we simply draw from the prior distribution $$p(\vec{w}) = \mathcal{N}(\vec{0}, I)$$ (where $$I$$ is the identity matrix), which we choose to be a multi-variate standard Gaussian. Then we gradually introduce new data points, compute the posterior and plot the distribution and the sampled models.

Wir können sehen, dass wir ohne Beobachtung der Datenpunkte im Durchschnitt eine Null vorhersagen, jedoch mit einer großen Varianz. Wenn wir mehr und mehr Datenpunkte beobachten, konvergiert die mittlere Vorhersage und die Standardabweichung nimmt ab. Intuitiv macht dies durchaus Sinn, da wir mit zunehmender Anzahl der beobachteten Daten immer sicherer in unseren Vorhersagen werden.

Lassen Sie uns nun etwas über die Mathematik hinter BLR herausfinden.

Der Satz von Bayes

Wie der Name BLR bereits andeutet, werden wir uns des Satzes von Bayes bedienen. Insbesondere werden wir die Formulierung des Satzes von Bayes für Wahrscheinlichkeitsdichten verwenden. Die Wahrscheinlichkeitsdichten $$p_X$$, $$p_{X|Y}$$ und $$p_{XY}$$ sind so definiert, dass

$$\begin{equation} \begin{aligned} \int_a^b p_X(x) dx &:= P( a < X < b ) \\ \int_a^b p_{X | Y} (x | y) dx &:= P( a < X < b | Y = y) \\ \int_c^d \int_a^b p_{XY} (x, y) dx dy &:= P( a < X < b , c < Y < d) \end{aligned} \end{equation}$$

für $$a \leq b$$ und $$c \leq d$$. Der Satz von Bayes für Wahrscheinlichkeitsdichten besagt, dass

$$\begin{equation} p_{Y|X}(y | x) = \frac{p_{X|Y}(x|y) p_Y(y)}{p_X(x)}. \label{eqBayesTheorem} \end{equation}$$

In gewissem Sinne erlaubt uns der Satz von Bayes, die bedingten Wahrscheinlichkeiten "umzudrehen"; alles, was wir brauchen, sind die unbedingten Verteilungen. Im Bayes'schen Rahmen werden die unbedingten Verteilungen als Priorverteilungen oder kurz Prioren bezeichnet, während die bedingten Verteilungen als Posteriorverteilungen bezeichnet werden. In der Praxis sieht der "Bayes'sche" Arbeitsablauf wie folgt aus:

  • Wir möchten $$p_{Y|X}(y|x)$$ berechnen, haben aber nur Zugriff auf $$p_{X|Y}(x|y)$$.

  • Daher müssen wir einen Prior $$p_Y(y)$$ erstellen, der unsere frühere "Überzeugung" über $$y$$ repräsentieren sollte.

  • Jetzt können wir das Produkt $$p_{X|Y}(x|y) p_Y(y)$$ berechnen, das für ein gegebenes $$x$$ eine (hoffentlich integrierbare) Funktion von $$y$$ ist.

  • Um eine integrierbare Funktion von $$y$$ in eine Wahrscheinlichkeitsdichte umzuwandeln, muss sie über eine Konstante $$N$$ (die von $$x$$ abhängen könnte) so normiert werden, dass $$\int p_{Y|X}(y|x) dy = \frac{1}{N} \int p_{X|Y}(x|y) p_Y(y) dy = 1$$, wobei die Integration über alle möglichen $$y$$ durchgeführt wird.

Der letzte Punkt impliziert

$$\begin{equation} N = \int p_{X|Y}(x|y) p_Y(y) dy = \int p_{XY}(x,y)dy = p_X(x), \label{eqNormalization} \end{equation}$$

was der Nenner des Bayes'schen Satzes ist. Kurz gesagt, wenn wir $$p_{Y|X}(y|x)$$ berechnen wollen, aber nur Zugriff auf $$p_{X|Y}(x|y)$$ haben, würden wir oft einfach ein vorheriges $$p_Y(y)$$ definieren und $$p_{Y|X}(y|x)$$ durch Normalisierung des Zählers des Bayes'schen Satzes erhalten. Die Berechnung der Normalisierung kann manchmal recht kostspielig sein, so dass uns diese Faulheit auf die Füße fallen könnte - es sei denn, es werden Schritte unternommen, um die Widerspenstigkeit des Normalisierungsintegrals in der obigen Formel zu mildern (z.B. Variationsinferenz).

Lassen Sie uns als nächstes die nicht-bayesianische lineare Regression genauer betrachten und erörtern, wie sie mit dem Bayes'schen Gegenstück zusammenhängt.  

Lineare Regression

Maximum Likelihood Estimator

Wenn Sie jemals ein kleines (oder manchmal sogar ein großes) Regressionsproblem gelöst haben, haben Sie höchstwahrscheinlich einen Maximum Likelihood Estimator (MLE) verwendet. Die Likelihood der Modellparameter ist definiert als

$$\begin{equation} \begin{aligned} L(\vec{w}) &= p(D|\vec{w}) \\ &= \prod_{i=1}^{N_D} p(y_i | \vec{x}_i, \vec{w}) p(\vec{x}_i) \label{eqLikelihood}. \end{aligned} \end{equation}$$

Da das Rauschen additiv ist, erhalten wir bequemerweise eine bedingte Verteilung

$$\begin{equation} p(y_i | \vec{x}_i, \vec{w}) = \mathcal{N}(\vec{w}^T \vec{x}_i, \sigma_\epsilon^2) \label{eqCondDistModel} \end{equation}.$$

Unser Ziel ist es nun, die Modellparameter $$\vec{w}_{MLE}^{\ast}$$ zu finden, die die oben definierte Wahrscheinlichkeit maximieren. Der Einfachheit halber verwenden wir oft die logarithmische Wahrscheinlichkeit $$LL$$. Da der Logarithmus strikt ansteigend ist, ändert die Verwendung des Logarithmus der Wahrscheinlichkeit die Lage des Maximums $$\vec{w}_{MLE}^{\ast}$$ nicht. Dies ergibt

$$\begin{equation} \begin{aligned} LL(\vec{w}) &= \log L(\vec{w}) \\ &= \log \prod_{i=1}^{N_D} p(y_i | \vec{x}_i, \vec{w}) p(\vec{x}_i) \\ &= \sum_{i=1}^{ND} \left[ \log p(y_i | \vec{x}_i, \vec{w}) + \log p(\vec{x}_i) \right] \\ &= \sum_{i=1}^{N_D} \left[ - \frac{(\vec{w}^T \vec{x}_i - y_i)^2}{2\sigma^2_\epsilon} - \log(\sqrt{2\pi \sigma^2_\epsilon} ) + \log p(\vec{x}_i) \right]. \end{aligned} \end{equation}$$

Der letzte Schritt folgt, da $$p(y_i | \vec{x}_i, \vec{w})$$ die Dichte einer Normalverteilung ist. Die Terme mit dem Logarithmus in der Summierung hängen nicht von den Modellparametern $$\vec{w}$$ ab, daher können sie bei der Optimierung ignoriert werden. Da auch die Skalierung des Optimierungsziels durch eine positive Zahl die Lage des Optimums nicht verändert, können wir außerdem den Faktor $$\frac{1}{2\sigma^2_\epsilon}$$ im ersten Term weglassen. Damit ergibt sich das Optimierungsproblem

$$\begin{equation} \begin{aligned} \vec{w}_{MLE}^{\ast} &= {\mathrm{argmax}}_{\vec{w}} L(\vec{w}) \\ &= {\mathrm{argmax}}_{\vec{w}} LL(\vec{w}) \\ &= {\mathrm{argmax}}_{\vec{w}} \sum_{i=1}^{N_D} -(\vec{w}^T \vec{x}_i - y_i)^2 \\ &= {\mathrm{argmin}}_{\vec{w}} \sum_{i=1}^{N_D} (\vec{w}^T \vec{x}_i - y_i)^2 . \end{aligned} \end{equation}$$

Das ist ein tolles Ergebnis. Der MLE für ein lineares Modell unter der Annahme eines additiven normalverteilten Rauschens erweist sich als Minimierung der Summe der quadrierten Fehler.

Maximum A-Posteriori

Der MLE ist schon ganz nett, manchmal könnte es jedoch vorteilhaft sein, vorherige Annahmen über die Modellparameter einzubeziehen, z.B. wenn wir nur eine kleine Menge verrauschter Daten haben. Dies kann mit Hilfe des Maximum A-Posteriori-Schätzers (MAP) erfolgen. Für den MLE haben wir versucht, den Parametersatz $$\vec{w}$$ so zu finden, dass die Wahrscheinlichkeit der beobachteten Daten $$p(D|\vec{w})$$ maximiert wird. Ein alternativer Ansatz bestünde darin, die Wahrscheinlichkeit $$p(\vec{w}|D)$$ der Parameter zu maximieren, nachdem die Daten $$D$$ beobachtet wurden. Da $$p(\vec{w}|D)$$ die posteriore Verteilung über die Parameter ist, hat die MAP ihren Namen von hier, d.h. wir wollen die Modellparameter $$\vec{w}^{\ast}_{MAP}$$ finden, die die posteriore Verteilung maximieren, wenn die von uns beobachteten Daten $$D$$ gegeben sind. Aber wie sieht $$p(\vec{w}|D)$$ aus?

Betrachten wir das Bayes'sche Theorem noch einmal:

$$\begin{equation} \begin{aligned} p(\vec{w} | D ) &= \frac{p(D|\vec{w})p(\vec{w})}{p(D)} \\ &= \frac{p(\vec{w}) \prod_{i=1}^{N_D} p(y_i | \vec{x}_i, \vec{w}) p(\vec{x}_i)}{\prod_{i=1}^{ND} p(y_i | \vec{x}_i) p(\vec{x}_i)} \\ &= \frac{p(\vec{w}) \prod_{i=1}^{ND} p(y_i | \vec{x}_i, \vec{w})}{\prod_{i=1}^{ND} p(y_i | \vec{x}_i)}. \end{aligned} \label{eqBayesForRegression} \end{equation}$$

Das Optimierungsproblem, das wir jetzt lösen müssen, kann wie folgt geschrieben werden:

$$\begin{equation} \begin{aligned} \vec{w}^{\ast}_{MAP} &={\mathrm{argmax}}_{\vec{w}} p(\vec{w} | D ) \\ &=  {\mathrm{argmax}}_{\vec{w}} \frac{p(\vec{w}) \prod_{i=1}^{N_D} p(y_i | \vec{x}_i, \vec{w})}{\prod_{i=1}^{N_D} p(y_i | \vec{x}_i)}. \end{aligned} \end{equation}$$

Da die $$p(y_i|\vec{x}_i)$$ nicht von den Modellparametern $$\vec{w}$$ abhängen, enden sie als einfache Faktoren, die die Lage von $$\vec{w}^{\ast}_{MAP}$$ nicht beeinflussen, deshalb können wir sie fallen lassen. Darüber hinaus können wir wieder die Monotonie des Logarithmus verwenden, um zu schreiben:

$$\begin{equation} \begin{aligned} \vec{w}^{\ast}_{MAP} &= {\mathrm{argmax}}_{\vec{w}} p(\vec{w}) \prod_{i=1}^{ND} p(y_i|\vec{x}_i, \vec{w}) \\ &= {\mathrm{argmax}}_{\vec{w}} \log \left[ p(\vec{w}) \prod_{i=1}^{N_D} p(y_i|\vec{x}_i, \vec{w}) \right] \\ &= {\mathrm{argmax}}_{\vec{w}} \left[ \log p(\vec{w}) + \sum_{i=1}^{ND} \log p(y_i|\vec{x}_i, \vec{w}) \right]. \end{aligned} \label{eqMAPObjective} \end{equation}$$

Wie wir sehen können, ist die Summierung über die Datenpunkte die gleiche wie im MLE, wir haben lediglich einen zusätzlichen Term $$\log p(\vec{w})$$ eingeführt, der es uns erlaubt, frühere Überzeugungen über die Modellparameter zu formulieren.

Für ein kurzes Beispiel nehmen wir an, dass unsere Parameter normal und unabhängig um den Ursprung mit einer Varianz $$\sigma^2_w$$ verteilt sind, d.h.

$$\begin{equation} \begin{aligned} p(\vec{w}) &= p(w_0)p(w_1) \\ &= \frac{1}{2\pi \sigma_w^2} \exp \left[ - \frac{w^2_0 + w^21}{2\sigma^2_w} \right] \\ &= \frac{1}{2\pi \sigma_w^2} \exp \left[ - \frac{\vec{w}^T \vec{w}}{2\sigma^2_w} \right]. \end{aligned} \label{eqMAPPrior} \end{equation}$$

Wenn wir die letzten beiden Gleichungen mit der oben abgeleiteten bedingten Verteilung kombinieren, erhalten wir

$$\begin{equation} \begin{aligned} \vec{w}^{\ast}_{MAP} &= {\mathrm{argmax}}_{\vec{w}} \left[ -\frac{\vec{w}^T \vec{w}}{2\sigma_w^2} - \sum_{i=1}^{N_D} \frac{(\vec{w}^T \vec{x}_i - y_i)^2}{2\sigma_\epsilon^2} \right], \end{aligned} \end{equation}$$

wo wir wieder parameterunabhängige Begriffe streichen. Da wir das Optimierungsziel mit jedem gewünschten positiven skalaren Wert skalieren können und grundsätzlich frei sind, $$\sigma_w^2$$ für unsere früheren Überzeugungen zu wählen, können wir das Optimierungsziel wie folgt umschreiben:

$$\begin{equation} \begin{aligned} \vec{w}^{\ast}_{MAP} &= {\mathrm{argmax}}_{\vec{w}} \left[ -\lambda \vec{w}^T \vec{w} - \sum_{i=1}^{N_D} (\vec{w}^T \vec{x}_i - y_i)^2 \right] \\ &= {\mathrm{argmin}}_{\vec{w}} \left[ \sum_{i=1}^{N_D} (\vec{w}^T \vec{x}_i - y_i)^2 + \lambda \vec{w}^T \vec{w}\right]. \end{aligned} \end{equation}$$

Wir sehen, dass MAP fast dasselbe Optimierungsproblem wie MLE mit einem zusätzlichen zusätzlichen Term löst, der durch irgendeinen Parameter $$\lambda$$ skaliert wird, für den wir einen vernünftigen Wert finden müssen. $$\lambda$$ quantifiziert im Wesentlichen, wie stark wir glauben, dass die Modellparameter nahe bei Null liegen. Leser mit einigen Kenntnissen im maschinellen Lernen werden erkennen, dass MAP dasselbe ist wie MLE mit L2-Regulierung.

Bayes'sche Lineare Regression

Wir haben gerade das Bayes'sche Theorem verwendet, um die Schätzung der Modellparameter mit einer Regularisierung zu rechtfertigen, aber wir verwenden immer noch Punktschätzungen der Modellparameter $$\vec{w}$$. Der Satz von Bayes könnte uns theoretisch nicht nur Zugriff auf das Maximum der posterioren Verteilung wie bei MAP geben, sondern uns auch erlauben, die Verteilung $$p(\vec{w} | D)$$ selbst zu verwenden. Dies wird uns erlauben, die Standardabweichung unserer Vorhersage zu schätzen. Wir werden konjugierte Prioren verwenden, um $$p(\vec{w} | D)$$ zu berechnen. Konjugierte Prioren sind Prioren $$p(\vec{w})$$, die zur gleichen Familie von Verteilungen gehören wie die späteren $$p(\vec{w} | D)$$. Wie in den obigen Beispielen gehen wir davon aus, dass es sich bei der Wahrscheinlichkeit $$p(y_i | \vec{x}_i, \vec{w})$$ und den früheren $$p(\vec{w})$$ um Normalverteilungen handelt. Eine Normalverteilung mal ein Normalverteilung ist immer noch eine Normalverteilung. Wie aus dem Satz von Bayes, der Formel für die Normalisierung und dem Satz von Bayes hervorgeht, ist der Prior der Trainingsdaten $$p(D)$$ im Wesentlichen eine Normalisierungskonstante und hat keinen Einfluss auf die allgemeine Form von $$p(\vec{w} | D)$$. Wir wissen also bereits, dass $$p(\vec{w} | D)$$ eine Normalverteilung sein wird. Für eine Normalverteilung müssen wir nur den das Mittel $$\vec{\mu}_w$$ und die Kovarianzmatrix $$\Sigma_w$$ von $$p(\vec{w} | D)$$ herausfinden und können dann daraus die Normalisierung ableiten.

Das war's fürs Erste. Im zweiten Teil werde ich die Details der Mathematik hinter der Bayes'schen linearen Regression erläutern.