基于R语言为同趋势时序数据拟合通用方程的技术请求

给定时序数据呈现初始峰值→平台波动→二次抬升的变化趋势，以下是在R中拟合可适配同类数据（趋势相似、数值水平不同）的方程方法：

示例数据 #

data <- structure(list(day = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253), Count = c(10, 50, 500, 425, 300, 400, 275, 98, 115, 79, 87, 114, 69, 105, 81, 82, 117, 87, 123, 81, 119, 97, 84, 124, 122, 53, 114, 95, 49, 95, 101, 114, 74, 120, 72, 61, 79, 59, 96, 95, 105, 53, 110, 69, 69, 79, 106, 52, 50, 98, 102, 107, 122, 108, 47, 68, 51, 114, 96, 102, 121, 113, 130, 134, 143, 144, 141, 139, 140, 142, 141, 125, 134, 130, 137, 139, 123, 138, 108, 133, 97, 122, 120, 110, 144, 121, 103, 127, 103, 100, 139, 138, 103, 105, 114, 142, 128, 141, 141, 122, 110, 125, 112, 98, 130, 116, 138, 120, 135, 143, 136, 145, 101, 120, 131, 119, 131, 116, 114, 143, 126, 102, 116, 106, 133, 110, 102, 141, 141, 132, 110, 95, 130, 133, 131, 128, 103, 111, 120, 140, 107, 114, 95, 113, 116, 131, 145, 144, 121, 111, 100, 145, 96, 130, 95, 119, 135, 127, 113, 105, 110, 102, 105, 116, 145, 115, 102, 120, 143, 140, 141, 132, 143, 136, 108, 106, 127, 112, 122, 118, 112, 96, 116, 141, 162, 168, 198, 156, 165, 180, 179, 166, 194, 194, 162, 199, 156, 193, 200, 160, 160, 187, 150, 185, 161, 183, 166, 167, 199, 159, 146, 195, 151, 161, 161, 162, 167, 193, 191, 181, 148, 200, 182, 164, 147, 182, 165, 165, 159, 163, 188, 154, 192, 157, 149, 163, 170, 151, 185, 168, 154, 164, 191, 169, 186, 157, 182, 195, 150, 145, 152, 188, 176)), row.names = c(NA, -253L), class = c("tbl_df", "tbl", "data.frame"))

拟合效果示例 #

拟合步骤 #

1. 加载依赖包 #

install.packages(c("nls2", "ggplot2")) # 首次安装需运行
library(nls2)
library(ggplot2)

2. 定义分段非线性模型 #

针对数据的三段趋势，定义连续分段模型（避免断点，适配不同数值水平）：

model <- Count ~ ifelse(day <= t1, 
                        peak * exp(-(day - peak_day)^2/(2*sigma^2)), # 初始高斯峰值段
                        ifelse(day <= t2, 
                               plateau, # 平台波动段（拟合时自动包含误差项）
                               plateau + slope*(day - t2))) # 二次抬升段

模型参数说明：

• peak：初始峰值高度
• peak_day：初始峰值对应的天数
• sigma：初始峰值的宽度
• t1：峰值段结束天数
• plateau：平台期均值
• t2：平台期结束天数
• slope：抬升段的斜率

3. 拟合模型 #

通过暴力搜索确定初始值范围，再用非线性最小二乘拟合：

# 设置初始值范围（可根据数据大致范围调整）
start_params <- list(peak = c(400, 600),
                     peak_day = c(2, 4),
                     sigma = c(1, 3),
                     t1 = c(5, 10),
                     plateau = c(80, 150),
                     t2 = c(150, 170),
                     slope = c(0.5, 2))

# 拟合模型
fit_result <- nls2(model, data = data, start = start_params, algorithm = "brute-force")

# 查看拟合参数
summary(fit_result)

4. 适配同类新数据 #

对于趋势相似但数值不同的同类数据，只需复用模型公式，调整初始值范围后重新拟合：

# 假设new_data为同类新时序数据
new_start_params <- list(peak = c(300, 700), # 根据新数据峰值范围调整
                         peak_day = c(1, 5),
                         sigma = c(1, 4),
                         t1 = c(4, 12),
                         plateau = c(70, 160), # 根据新数据平台范围调整
                         t2 = c(140, 180),
                         slope = c(0.3, 2.5))

new_fit <- nls2(model, data = new_data, start = new_start_params, algorithm = "brute-force")

5. 可视化拟合结果 #

# 生成拟合值
data$fitted_count <- predict(fit_result, data)

# 绘制拟合图
ggplot(data, aes(x = day)) +
  geom_point(aes(y = Count), alpha = 0.5, size = 1) +
  geom_line(aes(y = fitted_count), color = "red", linewidth = 1) +
  labs(x = "天数", y = "计数", title = "时序数据拟合效果") +
  theme_minimal()

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