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、古いバージョン）もあります。より広義には、packagemgcvは同じことを行いますが、異なるスプライン ベースのアプローチを使用します。

あなたがやっていることについては、グリッド上の予測に有限差分gam::gamまたは有限差分を使用します。mgcv::gam前者のみが実際のローカル回帰に基づいていますが、どちらも尋ねられている質問に答えています。

車輪を再発明する必要はないと思います。さらに重要なことは、既存のパッケージを使用するということは、エンドポイントでのバイアスや曲線の転換点でのバイアスなどの問題を考慮に入れることになるということです (ローカルの線形フィットは、ローカルの最大値/最小値に偏ります)。使用する重み付けスキーム。クロスバリデーションなど、モデルの構築とチェックのための標準ツールを利用することもできます。

于 2013-06-26T07:36:24.240 に答える
1

lmアイデアは、応答のマトリックスで 1 回だけ呼び出すことです。これは 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 に答える