Search code examples
rstatistics

Plotting the harmonic as a perturbation of the mean

I'm currently trying to plot the first and second harmonic (as depicted in Ramsay - FDA with R and MATLAB). Ideally, it should look like the image bellow - enter image description here So far, this is the set up I've got -

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)

For the mean fda script and the joint 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)

Could anyone point out how I add the variation of the harmonics to the mean? Or even a better way to re-write the code I have now. Thanks in advance!

Solution

  • So apparently, there is a predefined command that already does the job -

    plot.pca.fd(handPCA, harm = 1, pointplot = FALSE, cycle = TRUE)

    Although from the documentation, it isn't exactly clear that the cycle parameter will achieve this result.