在R中对年均温度数据进行正弦曲线拟合的问题求助

问题背景 #

已有45年每日温度CSV数据，已转换为年度维度（年度最高温、最低温、平均温），需用正弦函数 y = A*sin(omega * t + phi) + C 对年度平均温度进行最优拟合，原尝试代码报错。

原代码错误分析 #

1. 索引错误：R向量索引从1开始，原循环i in 0:45调用TMAX[i]会取空值，导致max()/min()返回Inf/-Inf
2. 数据不匹配：nls调用时传入原始每日数据data，但拟合目标是年度平均温度TMEAN，数据维度不对应
3. 初始参数不合理：omega直接设为2*pi未结合年度周期逻辑，且初始参数组合导致奇异梯度矩阵，无法进行拟合迭代

修正方案与代码实现 #

步骤1：正确计算初始参数 #

无需循环，直接用年度温度的全局最值计算初始A（振幅）和C（垂直偏移）：

# 计算初始A和C
A_init <- (max(TMAX) - min(TMIN)) / 2
C_init <- (max(TMAX) + min(TMIN)) / 2

步骤2：构造拟合专用数据集 #

将年度时间序列t和年度平均温度TMEAN整理为数据框，确保nls能正确识别变量：

fit_data <- data.frame(t = t, y = TMEAN)

步骤3：稳定拟合（两种可选方式） #

方式1：优化初始参数后用基础nls

调整omega为年度周期对应的值（每年1个周期，即2*pi/1），设置合理的phi初始值：

# 设置初始参数
start_params <- list(A = A_init, omega = 2*pi, phi = 0, C = C_init)

# 执行拟合
res_nls <- nls(y ~ A*sin(omega*t + phi) + C, 
               data = fit_data, 
               start = start_params,
               control = list(maxiter = 5000))

方式2：使用minpack.lm包的nlsLM（更稳定，避免奇异梯度）

如果基础nls仍报错，推荐用minpack.lm包的nlsLM，它对初始参数鲁棒性更强：

# 安装并加载包
install.packages("minpack.lm")
library(minpack.lm)

# 拟合
res_lm <- nlsLM(y ~ A*sin(omega*t + phi) + C, 
                data = fit_data, 
                start = start_params,
                control = list(maxiter = 5000))

步骤4：查看拟合结果与可视化 #

# 查看拟合参数
summary(res_nls) # 或summary(res_lm)

# 生成拟合曲线数据
fit_vals <- predict(res_nls, newdata = fit_data)

# 绘制原始数据与拟合曲线
plot(t, TMEAN, col = "purple", pch = 16, main = "年度平均温度与正弦拟合曲线",
     xlab = "年份", ylab = "平均温度(°C)")
lines(t, fit_vals, col = "red", lwd = 2)
legend("topright", legend = c("原始数据", "拟合曲线"), col = c("purple", "red"),
       pch = c(16, NA), lwd = c(NA, 2))

完整可运行代码 #

# 读取数据（替换为你的文件路径）
data <- read.csv("your_temperature_data.csv")

# 转换为年度维度
L <- nrow(data)
N <- length(seq(1, L, by = 365))
t <- seq(0, N-1, 1)
TMAX <- as.numeric(data[seq(1, L, by = 365), 2])
TMIN <- as.numeric(data[seq(1, L, by = 365), 3])
TMEAN <- as.numeric(data[seq(1, L, by = 365), 4])

# 计算初始参数
A_init <- (max(TMAX) - min(TMIN)) / 2
C_init <- (max(TMAX) + min(TMIN)) / 2
start_params <- list(A = A_init, omega = 2*pi, phi = 0, C = C_init)

# 构造拟合数据
fit_data <- data.frame(t = t, y = TMEAN)

# 执行拟合（推荐用nlsLM）
install.packages("minpack.lm")
library(minpack.lm)
res_lm <- nlsLM(y ~ A*sin(omega*t + phi) + C, 
                data = fit_data, 
                start = start_params,
                control = list(maxiter = 5000))

# 输出结果
summary(res_lm)

# 可视化
fit_vals <- predict(res_lm, newdata = fit_data)
plot(t, TMEAN, col = "purple", pch = 16, main = "Mesa, AZ 年度平均温度拟合",
     xlab = "年份", ylab = "平均温度(°C)")
lines(t, fit_vals, col = "red", lwd = 2)

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