回归中‘cumulated’期望值向量化及基于历史时点的y期望值计算

嘿，我明白你的需求了：要针对每个时间点，用该点之前所有历史数据拟合y ~ x1 + x2的回归模型，然后预测当前时间点y的期望值，还要实现高效的向量化处理对吧？

先确认下你的数据生成代码，没问题：

set.seed(42)
dat <- data.frame(t=1:1000,x1=runif(1000,1,10),x2=round(runif(1000,0,1)))
dat$y <- 8*dat$x1 - 5*dat$x2 + rnorm(1000)

下面给你两种高效的实现方案，其中第二种是真正的向量化/高性能处理：

方案1：基础R增量QR分解（高效循环） #

线性回归的核心是求解X'Xβ = X'y，我们可以用QR分解的增量更新来避免每次重新拟合整个模型，比显式循环调用lm快得多：

# 构造带截距项的设计矩阵
X <- cbind(1, dat$x1, dat$x2)
y <- dat$y

# 初始化存储预测值的向量，第一个时间点无历史数据，设为NA
y_hat <- numeric(nrow(dat))
y_hat[1] <- NA

# 从第一个观测开始初始化QR分解
qr_obj <- qr(X[1,,drop=FALSE])
beta <- solve.qr(qr_obj, y[1])

# 递推更新QR分解并计算预测值
for (i in 2:nrow(dat)) {
  # 把第i-1个观测加入到QR分解中（用前i-1个数据拟合）
  qr_obj <- qr.update(qr_obj, X[i-1,,drop=FALSE], y[i-1])
  # 求解当前回归系数
  beta <- solve.qr(qr_obj, qr.qty(qr_obj, y[1:(i-1)]))
  # 预测第i个时间点的y期望值
  y_hat[i] <- X[i,,drop=FALSE] %*% beta
}

# 合并到原数据中查看结果
dat$y_hat <- y_hat
head(dat)

方案2：用rollRegres包实现向量化扩展窗口回归 #

如果追求极致的效率（尤其是样本量更大时），推荐使用rollRegres包，它底层用C实现了增量QR分解，支持扩展窗口的滚动回归，完全避免了R层面的显式循环，属于真正的向量化处理：

首先安装包：

install.packages("rollRegres")

然后运行代码：

library(rollRegres)

# 执行扩展窗口滚动回归：
# - 拟合数据用前999个观测（对应预测第2-1000个时间点）
# - width参数设为1:999，表示预测第i个点时，用前i-1个数据拟合
# - newdata指定要预测的x1/x2（即第2-1000个时间点的特征）
fit <- roll_regres(
  y ~ x1 + x2,
  data = dat[-nrow(dat), ],
  width = 1:999,
  newdata = dat[-1, ]
)

# 把预测值合并到原数据，第一个时间点设为NA
dat$y_hat <- c(NA, fit$predictions)

# 查看前几个结果
head(dat)

为什么不推荐循环调用lm？ #

如果你尝试用for循环每次调用lm拟合前i-1个数据，效率会非常低——因为每次都要重新计算整个模型的矩阵分解，对于1000个样本来说会慢很多，上面两种方案都是通过增量更新来避免重复计算，效率提升明显。

内容的提问来源于stack exchange，提问作者bumblebee