5

任务：在移动窗口中找到最佳线性拟合的斜率（例如，最小化误差方差）。x 值是等距的，例如随着时间的推移自动测量。

问题：性能是一个问题，因为它需要对许多数据集重复。

天真的实现：循环 y 值。

#some data
x <- 0:(8*60)
set.seed(42)
y <- -x^2*0.01+x*20+rnorm(8*60+1,mean=300,sd=50)

plot(y~x,pch=".")

optWinLinFit0 <- function(x,y,win_length) {
  xfit <- x[seq_len(win_length)]
  xfit <- xfit-min(xfit)
  #regression on moving window
  res <- lapply(seq_len(length(x)-win_length),function(i,x,y) {
    y <- y[seq_len(win_length)+i-1]
    list(y=y,fit = lm.fit(cbind(1,xfit),y))    
  },x=x, y=y)
  #find fit with smallest sigma^2
  winner <- which.min(sapply(res,function(x) 1/(win_length-2)*sum(x$fit$residuals^2)))

  y <- res[[winner]]$y
  #return fit summary and predicted values
  list(n=winner,summary=summary(lm(y~xfit)),
       dat=data.frame(x=x[-seq_len(winner-1)][seq_len(win_length)],
                      y=y,
                      ypred=res[[winner]]$fit$fitted.values))
}
res0 <- optWinLinFit0(x,y,180)


lines(ypred~x,data=res0$dat,col="red",lwd=2)

红线给出了移动窗口位置的拟合值，其中误差方差最小： 在此处输入图像描述

任何想法如何更快地做到这一点？

4

2 回答 2

2

你基本上是在做一个内核回归。有很多为此设计的功能和包：KernSmooth，gam并locfit浮现在脑海中。在基础 R 中，也有loess（并且lowess，旧版本）。更广泛地说，包mgcv做同样的事情，但使用不同的基于样条的方法。

对于您正在做的事情，我会使用gam::gamormgcv::gam并在网格上的预测中使用有限差分。只有前者基于实际的局部回归，但它们都回答了所提出的问题。

我认为没有必要重新发明轮子。更重要的是，使用现有的包意味着您将考虑端点偏差和曲线转折点等问题（局部线性拟合将偏向于局部最大值/最小值）；使用的加权方案；等。您还可以利用标准工具进行模型构建和检查，例如交叉验证等。

于 2013-06-26T07:36:24.240 回答
1

这个想法是lm只用一个响应矩阵调用一次。这速度快了 2 倍，但假设 y 值不为零。如果可能有零值，您可以检查它并optWinLinFit0用作后备。

optWinLinFit1 <- function(x,y,win_length) {
  xfit <- x[seq_len(win_length)]
  xfit <- xfit-min(xfit)

  #get all windows of values in one matrix
  mat <- outer(y,rep(1,length(y)))

  require(Matrix)
  mat <- band(mat,k1=0,k2=win_length-1)
  mat <- as.matrix(mat)
  mat <- mat[,-(1:win_length-1)]
  nc <- ncol(mat)
  mat <- matrix(mat[mat!=0],ncol=nc)

  #regression with response matrix
  fit <- lm.fit(cbind(1,xfit),mat)

  #find fit with smallest sigma^2
  winner <- which.min(1/(win_length-2)*colSums(fit$residuals^2))

  y <- mat[,winner]
  #return fit summary and predicted values
  list(n=winner,
       summary=summary(lm(y~xfit)),
       dat=data.frame(x=x[-seq_len(winner-1)][seq_len(win_length)],
                      y=y,
                      ypred=fit$fitted.values[,winner])
  )
}

all.equal(res0$ypred,res1$ypred)
#[1] TRUE

library(microbenchmark)
microbenchmark(optWinLinFit0(x,y,180),optWinLinFit1(x,y,180),times=10)
# Unit: milliseconds
#                     expr      min       lq   median       uq      max neval
# optWinLinFit0(x, y, 180) 30.90678 31.73952 31.83930 35.61465 35.90352    10
# optWinLinFit1(x, y, 180) 12.76270 14.70842 15.70562 16.06347 17.41174    10
于 2013-06-25T16:25:55.767 回答