ر.ggplot2.إعادة إنشاء منحنى سلس من طريقة stat_smooth

Question

أحاول إعادة إنشاء الإجراء المستخدم لتقدير الوظيفة السلسة من stat_smooth في حزمة ggplot2.لنأخذ مثالا:

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.

الخط الناعم لا يناسب البيانات حقًا، لكنه لا يهم.الآن، دعونا نعيد إنشاء نفس الرسم البياني من البداية.وفق http://www.inside-r.org/r-doc/stats/loess نحتاج إلى استخدام نواة الترجيح ثلاثية التكعيبات ومتعددة الحدود من الدرجة 2 (افتراضيًا).لقد وجدت هذا http://www.maths.manchester.ac.uk/~peterf/MATH38011/NPR%20local%20Linear%20Estimator.pdf مقالة تصف كيفية تقدير وظيفة اللوس الملساء.أحاول إعادة إنشاء هذه الطريقة واستخدامها على بياناتي:

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
}

لست متأكدًا مما يجب أن أستخدمه للمعلمة ح, ، ولكن وفقًا لكتاب لدي "بالنسبة للنواة ثلاثية المكعبات ذات العرض المتري، ح هو نصف قطر منطقة الدعم." ومن ثم فإن أول شيء أحاوله هو:

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)

ومن الواضح أن هذه ليست نفس الوظيفة.لقد تم استخدام قيم مختلفة من ح لكن لم أتمكن من إعادة إنشاء نفس المنحنى، مما جعلني أعتقد أنني ارتكبت خطأً في مكان آخر.

Answer 1

ال stat_smooth(...) وظيفة في ggplot تقوم الحزمة فقط بتمرير بياناتك (التي يُحتمل أن تكون فرعية) إلى ملف loess(...) الوظيفة، كما يمكن توضيحه هنا:

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")

ال توثيق ل loess(...) يوفر مراجع لطريقة الحساب، على وجه التحديد، هنا.