R.gggrundstück2.Neuerstellung einer glatten Kurve aus der stat_smooth-Methode
Frage
Ich versuche, die Prozedur neu zu erstellen, die zum Schätzen der glatten Funktion aus stat_smooth im ggplot2-Paket verwendet wird.Nehmen wir ein Beispiel:
library(ggplot2)
n <- 100
X <- runif(n)*8
Y <- sin(3*X) + cos(X^2) + rnorm(n, 0, 0.5)
myData <- as.data.frame(cbind(X, Y))
p <- ggplot(myData, aes(y=Y, x=X)) +
stat_smooth(se = FALSE, size = 2) +
geom_point(size = 1)
p
geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
Die glatte Linie passt nicht wirklich zu den Daten, aber das spielt keine Rolle.Erstellen wir nun dasselbe Diagramm von Grund auf neu.GEM http://www.inside-r.org/r-doc/stats/loess wir müssen einen trikubischen Gewichtungskern und ein Polynom vom Grad 2 verwenden (standardmäßig).Ich habe das gefunden http://www.maths.manchester.ac.uk /~peterf/MATH38011/NPR%20lokal%20Linear%20Estimator.pdf artikel, der beschreibt, wie man die glatte Lössfunktion schätzt.Ich versuche, diese Methode neu zu erstellen und für meine Daten zu verwenden:
Dfct <- function(t){
if (abs(t) <= 1)
((1-abs(t)^3)^3) else
0
}
K_h <- function(x_0, x){
f_hat <- NULL
Dfct(abs(x - x_0)/h)
}
m_hat_loess <- function(X, Y){
e_1 <- c(1, 0, 0)
m_hat <- NULL
for(i in 1:length(X)){
K_h_vector <- NULL
for(j in 1:length(X)){
K_h_vector <- c(K_h_vector, K_h(X[i], X[j]))
}
X_0 <- cbind(rep(1, length(X)), (X - X[i]), (X - X[i])^2)
W <- diag(K_h_vector)
m_hat <- c(m_hat,
t(e_1)%*% solve(t(X_0)%*%W%*%X_0) %*% (t(X_0)%*%W%*%Y)
)
}
m_hat
}
Ich bin mir nicht sicher, was ich für Parameter verwenden soll h, aber laut einem Buch habe ich "Für Tri-Cube-Kernel mit metrischer Breite, h ist der Radius des Stützbereichs." Daher ist das erste, was ich versuche, ist:
h <- (max(X)-min(X))/2
Y_hat <- m_hat_loess(X, Y)
tempData <- as.data.frame(cbind(X, Y_hat))
ggplot(tempData , aes(x=X, y=Y_hat)) +
geom_line(size = 2)
Dies ist eindeutig nicht die gleiche Funktion.Ich habe verschiedene Werte von verwendet h konnte aber nicht die gleiche Kurve nachbilden, was mich glauben lässt, dass ich woanders einen Fehler gemacht habe.
Lösung
Der stat_smooth(...)
funktion in der ggplot
paket übergibt lediglich Ihre Daten (möglicherweise Teilmengen) an die loess(...)
funktion, wie hier demonstriert werden kann:
library(ggplot2)
set.seed(1)
n <- 100
X <- runif(n)*8
Y <- sin(3*X) + cos(X^2) + rnorm(n, 0, 0.5)
myData <- data.frame(X,Y)
fit <- loess(Y~X,data=myData)
myData$pred <- predict(fit)
ggplot(myData, aes(X,Y))+
geom_point()+
stat_smooth(se=F, size=3)+
geom_line(aes(X,pred),colour="yellow")
Der Dokumentation für loess(...)
enthält Verweise auf die Berechnungsmethode, insbesondere, hier.