Modelle zur Korrelations- und Regressionsanalyse
Moderator: EDi
Gwinzy
Beiträge: 8 Registriert: Do Okt 10, 2019 9:29 pm
Beitrag
von Gwinzy » Mi Jun 03, 2020 2:33 pm
Hallo,
ich soll eine kleinste-Quadrate Schätzung programmieren.
Meine Funktion ist die folgende:
Code: Alles auswählen
attach(swiss)
meinlogit <- function (x,y)
{
(solve(t(x)%*%x)[1,1]*t(x))%*%y
}
meinlogit(Education, Fertility)
Allerdings ist meine Ergebnis nicht das gleiche wie das von
Kann mir da vielleicht jemand helfen?
Liebe Grüße
Gwinzy
Zuletzt geändert von
jogo am Mi Jun 03, 2020 3:31 pm, insgesamt 1-mal geändert.
ruedi_br
Beiträge: 159 Registriert: Do Mär 01, 2018 3:53 pm
Beitrag
von ruedi_br » Mi Jun 03, 2020 2:55 pm
etwas output - der beiden Funktionen und z.B. "str" des Datensatzes - wäre als Input hilfreich.
VG
Ruedi
fortune(111)
jogo
Beiträge: 2085 Registriert: Fr Okt 07, 2016 8:25 am
Beitrag
von jogo » Mi Jun 03, 2020 3:34 pm
Hallo GwinZy,
Gwinzy hat geschrieben: ↑ Mi Jun 03, 2020 2:33 pm
Code: Alles auswählen
attach(swiss)
meinlogit <- function (x,y)
{
(solve(t(x)%*%x)[1,1]*t(x))%*%y
}
meinlogit(Education, Fertility)
Allerdings ist meine Ergebnis nicht das gleiche wie das von
Der Teil
[1,1] ist falsch. Richtig muss es heißen:
Code: Alles auswählen
meinlogit <- function (x,y)
{
(solve(t(x)%*%x)*t(x))%*%y
}
Dies kann noch verkürzt werden zu:
Code: Alles auswählen
meinlogit <- function (x,y) solve( t(x)%*%x, t(x)%*%y )
bzw.
Code: Alles auswählen
meinlogit <- function (x,y) solve( crossprod(x), crossprod(x, y))
Diese händische Berechnung ist von der numerischen Stabilität viel schlechter als der Algorithmus, der bei
lm() eingesetzt wird.
Gruß, Jörg
Gwinzy
Beiträge: 8 Registriert: Do Okt 10, 2019 9:29 pm
Beitrag
von Gwinzy » Mi Jun 03, 2020 4:16 pm
Hallo Ruedi und Jörg,
danke, aber ich habe jetzt eine andere Formel benutzt, mit der ich auf das richtige Ergebnis gekommen bin.
jogo
Beiträge: 2085 Registriert: Fr Okt 07, 2016 8:25 am
Beitrag
von jogo » Mi Jun 03, 2020 7:19 pm
Hallo Gwinzy,
die Formel interessiert mich.
Vielleicht kenne ich sie noch nicht.
Zeig doch mal bitte!
Gruß, Jörg
Gwinzy
Beiträge: 8 Registriert: Do Okt 10, 2019 9:29 pm
Beitrag
von Gwinzy » Fr Jun 05, 2020 4:03 pm
Code: Alles auswählen
attach(swiss)
meinlogit <- function (x,y)
{ idx <- (1:length(x))
nenner <- 0
zaehler <- 0
for(i in idx){
zaehler <- zaehler + (x[i]-mean(x))*(y[i]-mean(y))
nenner <- nenner + (x[i]-mean(x))^2
}
zaehler/nenner
}
/code]
jogo
Beiträge: 2085 Registriert: Fr Okt 07, 2016 8:25 am
Beitrag
von jogo » Fr Jun 05, 2020 9:39 pm
Hallo Gwinzy,
vielen Dank!
Kannst Du mir trotzdem zeigen, wie diese Funktion jetzt die gleichen Ergebnisse liefert wie
(wenigstens bezüglich der Koeffizienten der Regression).
Wie rufst Du diese Funktion im Zusammenhang mit den Daten auf?
Bei mir sieht es so aus mit meiner Variante:
Code: Alles auswählen
X <- cbind(1, swiss$Education)
Y <- swiss$Fertility
solve( crossprod(X), crossprod(X, Y))
lm(Fertility ~ Education, data=swiss)
Gruß, Jörg
Gwinzy
Beiträge: 8 Registriert: Do Okt 10, 2019 9:29 pm
Beitrag
von Gwinzy » Sa Jun 06, 2020 7:43 pm
Code: Alles auswählen
attach(swiss)
meinlogit <- function (x,y)
{ idx <- (1:length(x))
nenner <- 0
zaehler <- 0
beta <- 0
for(i in idx){
zaehler <- zaehler + (x[i]-mean(x))*(y[i]-mean(y))
nenner <- nenner + (x[i]-mean(x))^2
}
beta[2] <- zaehler/nenner
beta[1] <- mean(y)-beta[2]*mean(x)
return (list(beta[1], beta[2]))
}
meinlogit(Education, Fertility)
summary(lm(Fertility ~ Education))
meinlogit gibt dann das beta0 und beta1 (oder alpha und beta, je nach Definition) für den KQS zurück.
Die Werte stimmen mit den "Estimate" Werten von summary(lm()) überein:
Aber deine Lösung ist natürlich deutlich eleganter
jogo
Beiträge: 2085 Registriert: Fr Okt 07, 2016 8:25 am
Beitrag
von jogo » So Jun 07, 2020 10:03 am
Hallo Gwinzy,
klar - jetzt sehe ich es auch.
Du hast die beiden Normalgleichungen für die einfache lineare Regression für die Lösung genommen.
Soweit ist alles ok.
Allerdings machst Du nicht von vektorisiertem Rechnen Gebrauch. Es geht auch ohne for-Schleife:
Code: Alles auswählen
zaehler <- sum((x-mean(x))*(y-mean(y)))
nenner <- sum( (x-mean(x))^2)
... und
attach() ist Teufelszeug. Bitte nicht verwenden!
viewtopic.php?f=7&t=5
Gruß, Jörg
Gwinzy
Beiträge: 8 Registriert: Do Okt 10, 2019 9:29 pm
Beitrag
von Gwinzy » So Jun 07, 2020 2:05 pm
Danke für die Tipps, so schaut es echt gleich viel schöner aus