用2-3周移动平均值替换时间序列中的异常降雨量值

下面是一套完整的实现方案，从异常值识别到替换全流程：

1. 自动识别异常值 #

先通过IQR方法标记异常值（比手动设阈值更通用，适合批量数据）：

# 计算四分位距范围
q1 <- quantile(tbl$rainfall, 0.25)
q3 <- quantile(tbl$rainfall, 0.75)
iqr <- q3 - q1
upper_bound <- q3 + 1.5*iqr  # 降雨量不会为负，无需设置下限

# 标记异常行
tbl <- tbl %>%
  mutate(is_outlier = rainfall > upper_bound)

2. 计算不含异常值的移动平均值 #

重点：计算移动平均时要排除异常值，避免被异常值带偏。这里用3周中心窗口（当前周±1周），如果要2周的话调整窗口范围即可：

方法一：用slider包（灵活，支持任意窗口） #

先安装slider包（如果没装），然后计算：

# install.packages("slider")
library(slider)

tbl <- tbl %>%
  mutate(rolling_mean_3w = slide_dbl(
    row_number(),
    .f = function(i) {
      # 划定窗口范围：当前行前后1行
      window <- max(1, i-1):min(nrow(tbl), i+1)
      # 取窗口内非异常值的平均值
      mean(tbl$rainfall[window][!tbl$is_outlier[window]], na.rm = TRUE)
    }
  ))

方法二：手动用lag/lead计算（无需额外包） #

如果不想装新包，直接用dplyr的函数手动实现：

tbl <- tbl %>%
  mutate(
    # 收集当前、前1、后1周的降雨量，排除异常值
    window_vals = pmap(list(lag(rainfall), rainfall, lead(rainfall)), 
                       function(x,y,z) c(x,y,z)[!c(lag(is_outlier), is_outlier, lead(is_outlier))]),
    # 计算窗口平均值
    rolling_mean_3w = map_dbl(window_vals, mean, na.rm = TRUE)
  ) %>%
  select(-window_vals)  # 删除临时列

3. 替换异常值 #

把标记为异常的降雨量替换成计算好的移动平均值：

tbl <- tbl %>%
  mutate(rainfall_clean = ifelse(is_outlier, rolling_mean_3w, rainfall))

完整可运行代码 #

把所有步骤整合，直接复制就能跑：

library(tidyverse)
library(lubridate)
library(slider)  # 用手动方法可注释此行

# 原始数据
week <- c("10-03-2021", "10-10-2021", "10-17-2021", "10-24-2021", "10-31-2021", "11-6-2021", "11-13-2021", "11-20-2021", "11-27-2021")
rainfall <- c(11, 12, 13, 1010, 1000, 16, 16, 18, 15)

tbl <- tibble(week = mdy(week), 
              rainfall = rainfall)

# 1. 标记异常值
q1 <- quantile(tbl$rainfall, 0.25)
q3 <- quantile(tbl$rainfall, 0.75)
iqr <- q3 - q1
upper_bound <- q3 + 1.5*iqr

tbl <- tbl %>%
  mutate(is_outlier = rainfall > upper_bound)

# 2. 计算3周移动平均（排除异常值）
tbl <- tbl %>%
  mutate(rolling_mean_3w = slide_dbl(
    row_number(),
    .f = function(i) {
      window <- max(1, i-1):min(nrow(tbl), i+1)
      mean(tbl$rainfall[window][!tbl$is_outlier[window]], na.rm = TRUE)
    }
  ))

# 3. 替换异常值
tbl <- tbl %>%
  mutate(rainfall_clean = ifelse(is_outlier, rolling_mean_3w, rainfall))

# 查看结果
print(tbl)

运行后，示例中的1010会被替换成13（前后仅前一周的13是正常值），1000会被替换成16（前后仅后一周的16是正常值），完全符合合理替换的需求。

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