如何在R中选择合适的sin()项拟合含多周期的时间序列

我想用sin()函数拟合带有周期性（波峰和波谷）的时间序列，但目前只能通过猜测周期（如1个月、2个月、1年、2年等）的方式尝试。请问R中是否存在可以估计时间序列中多个周期的函数？

以下是我手动猜测周期后拟合的示例代码（下图红线为拟合结果），如何找到具有合适周期的sin()项？

t <- 1:365
y <- c(-1,-1.3,-1.6,-1.8,-2.1,-2.3,-2.5,-2.7,-2.9,-3,-2,-1.1,-0.3,0.5,1.1,1.6,2.1,2.5,2.8,3.1,3.4,3.7,4.2,4.6,5,5.3,5.7,5.9,6.2,5.8,5.4,5,4.6,4.2,3.9,3.6,3.4,3.1,2.9,2.8,2.6,2.5,2.3,1.9,1.5,1.1,0.8,0.5,0.2,0,-0.1,-0.3,-0.4,-0.5,-0.5,-0.6,-0.7,-0.8,-0.9,-0.8,-0.6,-0.3,-0.1,0.1,0.4,0.6,0.9,1.1,1.3,1.5,1.7,2.1,2.4,2.7,3,3.3,3.5,3.8,4.3,4.7,5.1,5.5,5.9,6.2,6.4,6.6,6.7,6.8,6.8,6.9,7,6.9,6.8,6.7,
       6.5,6.4,6.4,6.3,6.2,6,5.9,5.7,5.6,5.5,5.4,5.4,5.1,4.9,4.8,4.6,4.5,4.4,4.3,3.9,3.6,3.3,3,2.8,2.6,2.4,2.6,2.5,2.4,2.3,2.3,2.2,2.2,2.3,2.4,2.4,2.5,2.5,2.6,2.6,2.4,2.1,1.9,1.8,1.6,1.4,1.3,1,0.7,0.5,0.2,0,-0.2,-0.4,-0.2,-0.1,0.1,0.1,0.1,0.1,0.1,0.1,0,0,-0.1,-0.1,-0.2,-0.2,-0.3,-0.3,-0.4,-0.5,-0.5,-0.6,-0.7,-0.7,-0.8,-0.8,-0.8,-0.9,-0.9,-0.9,-1.3,-1.6,-1.9,-2.1,-2.3,-2.6,-2.9,-2.9,-2.9,-2.9,
       -2.9,-3,-3,-3,-2.8,-2.7,-2.5,-2.4,-2.3,-2.2,-2.1,-2,-2,-1.9,-1.9,-1.8,-1.8,-1.8,-1.9,-1.9,-2,-2.1,-2.2,-2.2,-2.3,-2.4,-2.5,-2.6,-2.7,-2.8,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.8,-2.8,-2.7,-2.7,-2.6,-2.6,-2.8,-3,-3.1,-3.3,-3.4,-3.5,-3.6,-3.5,-3.4,-3.3,-3.3,-3.2,-3,-2.9,-2.8,-2.8,-2.7,-2.6,-2.6,-2.6,-2.5,-2.6,-2.7,-2.8,-2.8,-2.9,-3,-3,-3,-3,-2.9,-2.9,-2.9,-2.9,-2.9,-2.8,
       -2.7,-2.6,-2.5,-2.4,-2.3,-2.3,-2.1,-1.9,-1.8,-1.7,-1.5,-1.4,-1.3,-1.5,-1.7,-1.8,-1.9,-2,-2.1,-2.2,-2.4,-2.5,-2.6,-2.7,-2.8,-2.8,-2.9,-3.1,-3.2,-3.3,-3.4,-3.5,-3.5,-3.6,-3.6,-3.5,-3.4,-3.3,-3.2,-3.1,-3,-2.7,-2.3,-2,-1.8,-1.5,-1.3,-1.1,-0.9,-0.7,-0.6,-0.5,-0.3,-0.2,-0.1,-0.3,-0.5,-0.6,-0.7,-0.8,-0.9,-1,-1.1,-1.1,-1.2,-1.2,-1.2,-1.2,-1.2,-0.8,-0.4,-0.1,0.2,0.5,0.8,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.6,0.3,0,-0.2,-0.5,-0.7,-0.8)

dt <- data.frame(t = t, y = y)
plot(x = dt$t, y = dt$y)

lm <- lm(y ~ sin(2*3.1416/365*t)+cos(2*3.1416/365*t)+
           sin(2*2*3.1416/365*t)+cos(2*2*3.1416/365*t)+
           sin(2*4*3.1416/365*t)+cos(2*4*3.1416/365*t)+
           sin(2*5*3.1416/365*t)+cos(2*5*3.1416/365*t)+
           sin(2*6*3.1416/365*t)+cos(2*6*3.1416/365*t)+
           sin(2*0.5*3.1416/365*t)+cos(2*0.5*3.1416/365*t),
         data = dt)

summary(lm)$adj.r.squared
plot(dt$y); lines(predict(lm), type = "l", col = "red")

解决方案 #

方法1：用傅里叶变换（FFT）识别显著周期 #

傅里叶变换能将时间序列转换到频域，定位能量最高的频率（对应周期），再用这些频率构建拟合模型，无需手动猜测。

# 计算傅里叶变换
fft_result <- fft(dt$y)
# 计算功率谱（振幅的平方，代表能量）
power <- Mod(fft_result)^2
# 生成频率序列，仅保留前半部分（频域结果对称）
freq <- seq_along(power)/length(power)
power_half <- power[2:floor(length(power)/2)]
freq_half <- freq[2:floor(length(power)/2)]

# 提取功率最高的前5个频率，转换为周期
top_freqs <- freq_half[order(power_half, decreasing = TRUE)[1:5]]
top_periods <- 1/top_freqs
cat("识别出的主要周期：", round(top_periods, 2), "\n")

# 用这些周期构建拟合公式
formula_terms <- lapply(top_freqs, function(f) {
  paste0("sin(2*pi*", f, "*t) + cos(2*pi*", f, "*t)")
})
formula_str <- paste("y ~", paste(unlist(formula_terms), collapse = " + "))
lm_fft <- lm(as.formula(formula_str), data = dt)

# 查看拟合效果
cat("调整后R平方：", round(summary(lm_fft)$adj.r.squared, 4), "\n")
plot(dt$y); lines(predict(lm_fft), type = "l", col = "blue")

方法2：使用forecast包自动生成傅里叶项 #

forecast包的fourier()函数可以快速生成指定数量的傅里叶项，结合线性模型拟合，大幅简化流程。

library(forecast)

# 生成10个傅里叶项（K值可根据拟合效果调整）
fourier_terms <- fourier(dt$t, K = 10)
dt_fourier <- cbind(dt, fourier_terms)

# 拟合模型
lm_fourier <- lm(y ~ ., data = dt_fourier)

# 查看结果
cat("调整后R平方：", round(summary(lm_fourier)$adj.r.squared, 4), "\n")
plot(dt$y); lines(predict(lm_fourier), type = "l", col = "green")

方法3：用非线性最小二乘（nls）直接估计周期参数 #

如果需要直接得到周期的具体估计值，可以用nls()构建非线性模型，让算法自动优化周期、振幅和相位参数。

# 定义多周期拟合模型
model_formula <- y ~ a0 + 
  a1*sin(2*pi*t/p1 + phi1) + 
  a2*sin(2*pi*t/p2 + phi2) +
  a3*sin(2*pi*t/p3 + phi3)

# 设置初始参数（只需大致范围，无需精确）
start_params <- list(a0 = mean(dt$y), 
                     a1 = 5, p1 = 365, phi1 = 0,
                     a2 = 2, p2 = 180, phi2 = 0,
                     a3 = 1, p3 = 90, phi3 = 0)

# 拟合模型
nls_fit <- nls(model_formula, data = dt, start = start_params)

# 查看估计的周期参数
summary(nls_fit)
# 绘制拟合曲线
plot(dt$y); lines(predict(nls_fit), type = "l", col = "purple")

方法对比 #

• 傅里叶变换：适合快速定位主要周期，数据量越大结果越可靠，适合初步筛选周期；
• forecast包傅里叶项：代码简洁，无需手动处理频域转换，适合快速建模；
• 非线性最小二乘：直接输出周期估计值，但需要合理的初始参数，否则可能不收敛。

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