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.

enter image description here

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)

enter image description here

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.

War es hilfreich?

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.

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top