Explizite Lösung für die lineare Regression mit zwei Prädiktoren

Ich habe einige Beispiele von Daten der Form und . Ich möchte eine Ebene an die Daten mit den kleinsten mittleren quadratischen Fehlern anpassen.$x,y$$z=f(x,y)$

$$z = Ax + By + C$$

Ich habe in Abschnitt 3 dieses Dokuments eine "Antwort" gefunden , die jedoch in Form einiger zu lösender Gleichungen verbleibt. Ich habe gerade die Fähigkeit, diese Gleichungen zu lösen, aber der Prozess wird so chaotisch, dass die Wahrscheinlichkeit, dass ich einen Fehler mache, ziemlich hoch ist. Sicherlich hat irgendwo jemand die vollständige Lösung in Langform geschrieben (sie könnte als "geschlossene Form" bezeichnet werden), und zwar in der Form

$$A=\ldots,$$ $$B=\ldots,$$ $$C=\ldots.$$

EDIT: Vielleicht ist "geschlossene Form" der falsche Ausdruck. Also lass mich klar sein. Ich möchte eine explizite Lösung für , und und keine Lösung, die mit "Wenn Sie diese Gleichungen lösen können, können Sie die Werte von , und " endet .$A$$B$$C$$A$$B$$C$

An anderer Stelle auf dieser Website finden Sie explizite Lösungen für die gewöhnliche Regression der kleinsten Quadrate

$$\mathbb{E}(z_i) = A x_i + B y_i + C$$

sind in Matrixform als erhältlich

$$(C,A,B)^\prime = (X^\prime X)^{-1} X^\prime z\tag{1}$$

wo $X$ ist die "Modellmatrix"

$$X = \pmatrix{1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ 1 & x_n & y_n}$$

und $z$ ist der Antwortvektor

$$z = (z_1, z_2, \ldots, z_n)^\prime.$$

Das ist eine vollkommen feine, explizite und berechenbare Antwort. Aber vielleicht gibt es ein zusätzliches Verständnis, das durch Inspektion der Koeffizienten herausgerissen werden kann. Dies kann erreicht werden, indem geeignete Einheiten ausgewählt werden, in denen die Variablen ausgedrückt werden sollen.

Die besten Einheiten für diesen Zweck zentrieren jede Variable auf ihren Mittelwert und verwenden ihre Standardabweichung als Maßeinheit. Lassen Sie die drei Mittel explizit sein$m_x, m_y,$ und $m_z$ und die drei Standardabweichungen sind $s_x, s_y,$ und $s_z$. (Es stellt sich heraus, dass es keine Rolle spielt, ob Sie durch teilen$n$ oder $n-1$bei der Berechnung der Standardabweichungen. Stellen Sie einfach sicher, dass Sie eine konsistente Konvention verwenden, wenn Sie einen zweiten Moment der Daten berechnen.) Die Werte der Variablen in diesen neuen Maßeinheiten sind

$$\xi_i = \frac{x_i - m_x}{s_x},\ \eta_i = \frac{y_i - m_y}{s_y},\ \zeta_i = \frac{z_i - m_z}{s_z}.$$

Dieser Prozess wird als Standardisierung der Daten bezeichnet. Die Variablen$\xi$, $\eta$, und $\zeta$sind die standardisierten Versionen der ursprünglichen Variablen$x$, $y$, und $z$.

Diese Beziehungen sind invertierbar:

$$x_i = s_x \xi_i + m_x,\ y_i = s_y \eta_i + m_y,\ z_i = s_z \zeta_i + m_z.$$

Einfügen dieser in die definierende Beziehung

$$\mathbb{E}(z_i) = C + Ax_i + By_i$$

und Vereinfachung der Erträge

$$\mathbb{E}(s_z \zeta_i + m_z) = C + A(s_x \xi_i + m_x) + B(s_y \eta_i + m_y).$$

Lösung für die Erwartung der abhängigen Variablen $\zeta_i$ ergibt

$$\mathbb{E}(\zeta_i) = \left(\frac{C + Am_x + Bm_y - m_z}{s_z}\right) + \left(\frac{A s_x}{s_z}\right) \xi_i + \left(\frac{B s_y}{s_z}\right) \eta_i.$$

Wenn wir diese Koeffizienten schreiben als $\beta_0, \beta_1, \beta_2$ jeweils können wir uns dann erholen $A, B, C$durch Vergleichen und Lösen. Für die Aufzeichnung gibt dies

$$A = \frac{s_z \beta_1}{s_x},\ B = \frac{s_z \beta_2}{s_y},\text{ and }C = s_z \beta_0 + m_z - A m_x - B m_y.\tag{2}$$

Der Punkt davon wird deutlich, wenn wir die neue Modellmatrix betrachten

$$\Xi = \pmatrix{1 & \xi_1 & \eta_i \\ 1 & \xi_2 & \eta_2 \\ \vdots & \vdots & \vdots \\ 1 & \xi_n & \eta_n}$$

und die neue Antwortmatrix $\zeta = (\zeta_1, \zeta_2, \ldots, \zeta_n)$, weil jetzt

$$\Xi^\prime \Xi = \pmatrix{n & 0 & 0 \\ 0 & n & n\rho \\ 0 & n\rho & n}$$

und

$$\Xi^\prime \zeta = (0, n\tau, n\upsilon)^\prime$$

wo $\rho$ ist der Korrelationskoeffizient $\frac{1}{n}\sum_{i=1}^n \xi_i \eta_i$, $\tau$ ist der Korrelationskoeffizient $\frac{1}{n}\sum_{i=1}^n \xi_i \zeta_i$, und $\upsilon$ ist der Korrelationskoeffizient $\frac{1}{n}\sum_{i=1}^n \eta_i \zeta_i$.

Um die normalen Gleichungen zu lösen $(1)$ wir können beide Seiten durch teilen $n$geben

$$\pmatrix{1 & 0 & 0 \\ 0 & 1 & \rho \\ 0 & \rho & 1}\pmatrix{\beta_0 \\ \beta_1 \\ \beta_2} = \pmatrix{0 \\ \tau \\ \upsilon} .$$

Was ursprünglich wie eine beeindruckende Matrixformel aussah, wurde auf einen wirklich elementaren Satz von drei simultanen Gleichungen reduziert. Unter der Voraussetzung$|\rho| \lt 1$ist seine Lösung leicht zu finden

$$\pmatrix{\hat\beta_0 \\ \hat\beta_1 \\ \hat\beta_2} = \frac{1}{1-\rho^2}\pmatrix{0 \\ \tau-\rho\upsilon \\ \upsilon-\rho\tau}.$$

Stecken Sie diese in die Koeffizienten in $(2)$ erzeugt die Schätzungen $\hat A, \hat B,$ und $\hat C$.

In der Tat wurde noch mehr erreicht:

• Es ist jetzt offensichtlich, warum die Fälle $|\rho|=1$ sind problematisch: Sie führen eine Division durch Null in die Lösung ein.

• Es ist ebenso offensichtlich, wie zu bestimmen ist, ob wann eine Lösung existiert $|\rho=1|$und wie man es erhält. Es wird existieren, wenn die zweite und dritte Normalgleichung in$\Xi$ sind redundant und werden einfach durch Ignorieren einer der Variablen erhalten $x$ und $y$ an erster Stelle.

• Wir können allgemein einen Einblick in die Lösung gewinnen. Zum Beispiel von$\hat\beta_0=0$In allen Fällen können wir den Schluss ziehen, dass die angepasste Ebene den Mittelwertpunkt durchlaufen muss$(m_x, m_y, m_z)$.

• Es ist nun offensichtlich, dass die Lösung in Bezug auf die ersten beiden Momente des trivariaten Datensatzes gefunden werden kann $(x, y, z)$. Dies wirft ein weiteres Licht auf die Tatsache, dass Koeffizientenschätzungen allein aus Mittelwerten und Kovarianzmatrizen ermittelt werden können .

• Weiterhin Gleichung $(2)$ zeigt, dass die Mittel nur zur Schätzung des Intercept-Terms benötigt werden $C$. Schätzungen der beiden Pisten$A$ und $B$ benötigen nur die zweiten Momente.

• Wenn die Regressoren nicht korreliert sind, $\rho=0$ und die Lösung besteht darin, dass der Achsenabschnitt Null ist und die Steigungen die Korrelationskoeffizienten zwischen der Antwort sind $z$ und die Regressoren $x$ und $y$ wenn wir die Daten standardisieren. Dies ist sowohl leicht zu merken als auch bietet einen Einblick in die Beziehung zwischen Regressionskoeffizienten und Korrelationskoeffizienten.

---

Wenn wir das alles zusammenfassen, finden wir das (außer in den entarteten Fällen) $|\rho|=1$) Die Schätzungen können geschrieben werden

$$\eqalign{ \hat A &= \frac{\tau - \rho\upsilon}{1-\rho^2} \frac{s_z}{s_x} \\ \hat B &= \frac{\upsilon - \rho\tau}{1-\rho^2} \frac{s_z}{s_y} \\ \hat C &= m_z -m_x \hat A - m_y \hat B. }$$

In diesen Formeln ist die $m_{*}$ sind die Beispielmittel, die $s_{*}$ sind die Standardabweichungen der Stichprobe und die griechischen Buchstaben $\rho, \tau,$ und $\upsilon$ repräsentieren die drei Korrelationskoeffizienten (zwischen $x$ und $y$, $x$ und $z$, und $y$ und $z$, beziehungsweise).

Bitte beachten Sie, dass diese Formeln nicht der beste Weg sind, um die Berechnungen durchzuführen. Sie alle beinhalten das Subtrahieren von Mengen, die von vergleichbarer Größe sein könnten, wie z$\tau-\rho\upsilon$, $\upsilon-\rho\tau$, und $m_z - (-m_x \hat A - m_y \hat B)$. Eine solche Subtraktion beinhaltet einen Genauigkeitsverlust. Die Matrixformulierung ermöglicht es numerischen Analysten, stabilere Lösungen zu erhalten, die so viel Präzision wie möglich bewahren. Aus diesem Grund haben Menschen selten Interesse an termingerechten Formeln. Der andere Grund, warum wenig Interesse besteht, ist, dass mit zunehmender Anzahl von Regressoren die Komplexität der Formeln exponentiell zunimmt und schnell zu unhandlich wird.

---

Als weiteren Beweis für die Richtigkeit dieser Formeln können wir ihre Antworten mit denen eines Standardlösers für kleinste Quadrate vergleichen, der lmFunktion in R.

#
# Generate trivariate data.
#
library(MASS)
set.seed(17)
n <- 20
mu <- 1:3
Sigma <- matrix(1, 3, 3)
Sigma[lower.tri(Sigma)] <- Sigma[upper.tri(Sigma)] <- c(.8, .5, .6)
xyz <- data.frame(mvrnorm(n, mu, Sigma))
names(xyz) <- c("x", "y", "z")
#
# Obtain the least squares coefficients.
#
beta.hat <- coef(lm(z ~ x + y, xyz))
#
# Compute the first two moments via `colMeans` and `cov`.
#
m <- colMeans(xyz)
sigma <- cov(xyz)
s <- sqrt(diag(sigma))
rho <- t(t(sigma/s)/s); rho <- as.vector(rho[lower.tri(rho)])
#
# Here are the least squares coefficient estimates in terms of the moments.
#
A.hat <- (rho[2] - rho[1]*rho[3]) / (1 - rho[1]^2) * s[3] / s[1]
B.hat <- (rho[3] - rho[1]*rho[2]) / (1 - rho[1]^2) * s[3] / s[2]
C.hat <- m[3] - m[1]*A.hat - m[2]*B.hat
#
# Compare the two solutions.
#
rbind(beta.hat, formulae=c(C.hat, A.hat, B.hat))

Die Ausgabe weist erwartungsgemäß zwei identische Schätzreihen auf:

         (Intercept)         x        y
beta.hat    1.522571 0.3013662 0.403636
formulae    1.522571 0.3013662 0.403636

Fand es...

$$\begin{align}A&=\frac{(\overline{zx}-\bar z\bar x)(\overline{y^2}-\bar y^2)-(\overline{zy}-\bar z\bar y)(\overline{xy}-\bar x\bar y)}{(\overline{x^2}-\bar x^2)(\overline{y^2}-\bar y^2)-(\overline{xy}-\bar x\bar y)^2}\\ B&=\frac{(\overline{x^2}-\bar x^2)(\overline{zy}-\bar z\bar y)-(\overline{zx}-\bar z\bar x)(\overline{xy}-\bar x\bar y)}{(\overline{x^2}-\bar x^2)(\overline{y^2}-\bar y^2)-(\overline{xy}-\bar x\bar y)^2}\\ C&=\bar z-A\bar x-B\bar y.\end{align}$$

Aber du hättest es verstehen sollen! Schauen Sie sich Ihre Entwurfsmatrix noch einmal an (erste Gleichung Seite 4). Elemente werden als (x, y) geschrieben.

Element (3,1) = Element (1,3)

Element (2,1) = Element (1,2)

Element (3,2) = Element (2,3)

und Element (3,3) ist m.

Ihre reguläre Entwurfsmatrix wird aufgelöst in (alle Beträge ersetzt)

$$\left( {\begin{array}{*{20}{c}} K&L&M\\ L&R&S\\ M&S&m \end{array}} \right)$$ Dies erleichtert das Lösen von Mach, insbesondere wenn Sie noch ein paar Schritte voraus denken. Falls Sie das getan haben, Bravo! Falls Sie es nicht getan haben, haben Sie nicht lange genug auf die Designmatrix gestarrt.