基于R语言FDA：如何将谐波变异叠加至均值并绘图？

我参考Ramsay所著《FDA with R and MATLAB》，尝试绘制一阶和二阶谐波，期望效果为包含均值曲线及叠加±1倍、±2倍一阶/二阶谐波变异的曲线集合（类似手部手写轨迹的变异展示图）。

已完成的初始代码如下：

平滑处理代码 #

handpartx <- handwrit[,,1]
handparty <- handwrit[,,2]
curv.Lfd <- int2Lfd(2)
lambda <- 1e-2
breaks <- seq(min(handwritTime), max(handwritTime), length.out = 100)
xbasis <- create.bspline.basis(range(handwritTime), norder = 6,
                               breaks = breaks)
ybasis <- create.bspline.basis(range(handwritTime), norder = 6,
                               breaks = breaks)

curvx.fdPar <- fdPar(xbasis, curv.Lfd, lambda)
curvy.fdPar <- fdPar(ybasis, curv.Lfd, lambda)

handxsmooth <- smooth.basis(handwritTime, handpartx, curvx.fdPar)
handysmooth <- smooth.basis(handwritTime, handparty, curvy.fdPar)
mean_x <- eval.fd(fdobj = mean.fd(handxsmooth$fd), evalarg = handwritTime)
mean_y <- eval.fd(fdobj = mean.fd(handysmooth$fd), evalarg = handwritTime)

联合PCA代码 #

zbasis <- smooth.basis(handwritTime, handwrit, curvx.fdPar)$fd
zbasis$fdnames[[1]] <- "Time"
zbasis$fdnames[[2]] <- "Replications"
zbasis$fdnames[[3]] <- list("X", "Y")

nharm <- 10
handPCA <- pca.fd(zbasis, nharm)

代码优化方案 #

首先可以简化重复代码，统一处理X/Y方向的平滑逻辑，避免冗余：

# 统一创建基础样条和正则化参数，X/Y共享相同配置
curv.Lfd <- int2Lfd(2)
lambda <- 1e-2
breaks <- seq(min(handwritTime), max(handwritTime), length.out = 100)
shared_basis <- create.bspline.basis(range(handwritTime), norder = 6, breaks = breaks)
shared_fdPar <- fdPar(shared_basis, curv.Lfd, lambda)

# 联合平滑X/Y二维数据，直接生成适配PCA的fd对象
hand_smooth <- smooth.basis(handwritTime, handwrit, shared_fdPar)
hand_fd <- hand_smooth$fd
# 标准化fd对象的名称，方便后续识别
hand_fd$fdnames <- list("Time", "Replications", c("X", "Y"))

# 计算均值函数
mean_fd <- mean.fd(hand_fd)

# 执行联合PCA
nharm <- 10
handPCA <- pca.fd(hand_fd, nharm)

叠加谐波变异至均值并绘制的实现 #

要实现目标效果，核心是提取主成分（谐波）并将其按标准差倍数叠加到均值曲线上，具体代码如下：

# 定义更密集的评估时间点，让曲线更平滑
eval_time <- seq(min(handwritTime), max(handwritTime), length.out = 200)

# 评估均值曲线的X/Y坐标
mean_x <- eval.fd(eval_time, mean_fd[, , 1])
mean_y <- eval.fd(eval_time, mean_fd[, , 2])

# 提取前2个主成分（一阶、二阶谐波）及对应标准差
top_harmonics <- handPCA$harmonics[, , 1:2]
# 主成分的标准差为方差的平方根（handPCA$values存储的是方差）
harm_sd <- sqrt(handPCA$values[1:2])

# 定义要叠加的变异倍数：±1、±2倍标准差
multipliers <- c(-2, -1, 1, 2)

# 初始化绘图，先绘制均值曲线
plot(mean_x, mean_y, type = "l", lwd = 2, col = "black", 
     xlab = "X坐标", ylab = "Y坐标", main = "均值曲线±谐波变异")

# 叠加一阶谐波的变异曲线
for (k in multipliers) {
  var_x <- mean_x + k * harm_sd[1] * eval.fd(eval_time, top_harmonics[, , 1])
  var_y <- mean_y + k * harm_sd[1] * eval.fd(eval_time, top_harmonics[, , 2])
  # 用不同颜色区分倍数，虚线表示一阶谐波变异
  lines(var_x, var_y, lty = 2, col = ifelse(abs(k)==2, "darkred", "darkblue"))
}

# 叠加二阶谐波的变异曲线
for (k in multipliers) {
  var_x <- mean_x + k * harm_sd[2] * eval.fd(eval_time, top_harmonics[, , 2])
  var_y <- mean_y + k * harm_sd[2] * eval.fd(eval_time, top_harmonics[, , 2])
  # 用不同颜色区分倍数，点虚线表示二阶谐波变异
  lines(var_x, var_y, lty = 3, col = ifelse(abs(k)==2, "darkgreen", "darkcyan"))
}

# 添加图例，清晰区分各曲线含义
legend("topright", legend = c("均值曲线", 
                              "一阶谐波±1σ", "一阶谐波±2σ",
                              "二阶谐波±1σ", "二阶谐波±2σ"),
       lty = c(1,2,2,3,3),
       col = c("black", "darkblue", "darkred", "darkcyan", "darkgreen"),
       bty = "n")

关键说明 #

• 联合处理X/Y数据：既减少重复代码，又保证平滑参数一致，和PCA输入格式匹配，降低出错概率
• 谐波变异幅度：pca.fd返回的values是各主成分的方差，取平方根得到标准差，以此控制变异的幅度，符合统计意义上的离散程度
• 叠加逻辑：均值曲线 + 倍数×标准差×谐波函数值，得到的是该主成分方向上的变异曲线，对应不同倍数的离散程度
• 可视化区分：用不同线型和颜色区分一阶/二阶谐波以及不同倍数的变异，让图形层次清晰，符合目标效果要求

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