R语言按SID分组计算Tor时序均值方差实现超限预测

现有包含100万行记录的数据框，需针对每个SID分组，沿时间维度计算指标Tor的均值与方差，以此预测Tor值超出限定范围（下限0.4、上限0.7）的时间节点。数据样例如下：

dat <- structure(list(timestamp = c("29-06-2021-06:00", "29-06-2021-06:01", 
"29-06-2021-06:02", "29-06-2021-06:03", "29-06-2021-06:04", "29-06-2021-06:05", 
"29-06-2021-06:06", "29-06-2021-06:07", "29-06-2021-06:08", "29-06-2021-06:09", 
"29-06-2021-06:10", "29-06-2021-06:11", "29-06-2021-06:12", "29-06-2021-06:13", 
"29-06-2021-06:14", "29-06-2021-06:15", "29-06-2021-06:16", "29-06-2021-06:17", 
"29-06-2021-06:18", "29-06-2021-06:19", "29-06-2021-06:20", "29-06-2021-06:21", 
"29-06-2021-06:22", "29-06-2021-06:23", "29-06-2021-06:24", "29-06-2021-06:25", 
"29-06-2021-06:26"), SID = c(301L, 351L, 304L, 357L, 358L, 302L, 
303L, 309L, 356L, 304L, 308L, 351L, 304L, 357L, 358L, 302L, 303L, 
352L, 307L, 353L, 304L, 308L, 352L, 307L, 304L, 354L, 356L), 
    Tor = c(0.70161919, 0.639416295, 0.288282073, 0.932362166, 
    0.368616626, 0.42175565, 0.409735918, 0.942170196, 0.381396521, 
    0.818102394, 0.659391671, 0.246387978, 0.196001777, 0.632630259, 
    0.66618385, 0.440625167, 0.639759498, 0.050001835, 0.775660271, 
    0.762934189, 0.516830196, 0.244674975, 0.38620466, 0.970792903, 
    0.752674581, 0.190366737, 0.56596405), Lowt = c(0L, 0L, 1L, 
    0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
    0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L), Hit = c(1L, 0L, 0L, 
    1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L)), class = "data.frame", row.names = c(NA, 
-27L))

head(dat)
#         timestamp SID       Tor Lowt Hit
#1 29-06-2021-06:00 301 0.7016192    0   1
#2 29-06-2021-06:01 351 0.6394163    0   0
#3 29-06-2021-06:02 304 0.2882821    1   0
#4 29-06-2021-06:03 357 0.9323622    0   1
#5 29-06-2021-06:04 358 0.3686166    1   0
#6 29-06-2021-06:05 302 0.4217556    0   0

• timestamp：样本记录时间
• SID：采集读数的部件ID，取值范围为301-310、351-360
• Tor：实际读数，数据类型为双精度浮点型<dbl>
• Lowt：二元变量，用于标识Tor读数低于下限
• Hit：二元变量，用于标识Tor读数低于上限

注：从样例数据看，Hit=1对应的Tor值均高于0.7上限，可核对字段标识逻辑是否存在笔误，避免后续统计出错。

时序场景下计算均值方差不能使用全量数据，否则会产生数据泄露，必须使用滚动/扩展窗口计算：即每个时间点的统计值仅能使用当前时间点及之前的历史数据，不能调用未来时间的读数。针对100万行的数据规模，优先选择高性能实现避免计算卡顿。

步骤1：数据预处理 #

先转换时间格式，按SID分组后沿时间升序排序，保证时序计算顺序正确：

library(data.table)
library(RcppRoll) # 基于Rcpp的高性能滚动计算函数，百万行数据秒级完成

# 转换为data.table高性能格式，解析时间字段
setDT(dat)
dat[, timestamp := as.POSIXct(timestamp, format = "%d-%m-%Y-%H:%M")]
# 每个SID组内按时间从早到晚排序
setorder(dat, SID, timestamp)

步骤2：分组计算时序均值与方差 #

根据业务需求选择窗口类型，两种常用方案如下：

# ---------- 方案1：固定大小滚动窗口 ----------
# 示例窗口大小为5，即每个时间点用前序4条+自身共5条历史数据计算统计值
# 窗口大小可根据采样频率调整，比如1分钟采样要预判10分钟风险可设为10-15
win_size <- 5
dat[, `:=`(
  roll_mean = RcppRoll::roll_mean(Tor, n = win_size, align = "right", fill = NA),
  roll_var = RcppRoll::roll_var(Tor, n = win_size, align = "right", fill = NA)
), by = SID]

# ---------- 方案2：无固定长度扩展窗口 ----------
# 每个时间点用该SID从首个观测点到当前点的全部历史数据计算统计值
dat[, `:=`(
  cum_mean = cummean(Tor),
  # 总体方差除以样本量n，样本方差除以n-1，按需调整
  cum_var = cumsum((Tor - cummean(Tor))^2)/seq_len(.N)
), by = SID]

步骤3：越限时间节点预测 #

采用工业场景常用的3σ原则判断：基于历史均值和方差计算置信区间，当置信区间碰触0.4/0.7阈值时，标记为即将越限，取每个SID首次触发预警的时间点即可：

# 以固定滚动窗口的统计值为例，计算95%置信区间（对应均值±2倍标准差，标准差=方差的平方根）
dat[, `:=`(
  pred_lower = roll_mean - 2*sqrt(roll_var),
  pred_upper = roll_mean + 2*sqrt(roll_var)
)]

# 提取每个SID首次触发越限预警的时间
first_warn <- dat[
  !is.na(roll_mean) & (pred_lower < 0.4 | pred_upper > 0.7),
  .(SID, first_warn_time = min(timestamp))
]

• 窗口大小不要过大或过小：窗口过大会对趋势变化反应迟钝，预警滞后；窗口过小容易被随机噪声干扰，产生误报，可结合历史越限数据回测选择最优窗口大小。
• 方差计算需区分总体方差和样本方差：如果窗口内数据是全部研究总体，计算时除以数据量n即可；如果是用窗口样本估计整体分布，除以n-1得到无偏估计更准确。
• 排序步骤不可省略：必须保证每个SID组内数据按时间升序排列，否则滚动计算的时序逻辑错误，结果完全无效。

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