R/AmpPhaseDecomp.R In fda: Functional Data Analysis

Documented in AmpPhaseDecomp

```AmpPhaseDecomp <- function(xfd, yfd, hfd, rng=xrng)
{
#  Computes the amplitude-phase decomposition for a registration.

#  Arguments:
#  XFD  ...  FD object for unregistered functions
#  YFD  ...  FD object for registered functions
#  HFD  ...  FD object for warping functions

#  Returns:
#  MS.amp ... mean square for amplitude variation
#  MS.pha ... mean square for amplitude variation
#  RSQR   ... squared correlation measure of prop. phase variation
#  C      ... constant C

xbasis  <- xfd\$basis
nxbasis <- xbasis\$nbasis
nfine   <- max(201,10*nxbasis)
xrng    <- xbasis\$rangeval

if (rng[1] < xrng[1] || rng[2] > xrng[2]) stop(
"RNG is not within the range of the other arguments.")
if (rng[1] >= rng[2]) stop("Elements of rng are not increasing.")
tfine   <- seq(rng[1],rng[2],len=nfine)
delta   <- tfine[2] - tfine[1]

Dhfine  <- eval.fd(tfine, hfd, 1)
xfine   <- eval.fd(tfine, xfd, 0)
yfine   <- eval.fd(tfine, yfd, 0)
mufine  <- apply(xfine, 1, mean)
etafine <- apply(yfine, 1, mean)

N       <- dim(xfine)[2]
rfine   <- yfine - outer(etafine,rep(1,N))

intetasqr <- delta*trapz(etafine^2)
intmusqr  <- delta*trapz(mufine^2)

covDhSy <- rep(0,nfine)
for (i in 1:nfine) {
Dhi        <- Dhfine[i,]
Syi        <- yfine[i,]^2
covDhSy[i] <- cov(Dhi, Syi)
}
intcovDhSy <- delta*trapz(covDhSy)

intysqr <- rep(0,N)
intrsqr <- rep(0,N)
for (i in 1:N) {
intysqr[i] <- delta*trapz(yfine[,i]^2)
intrsqr[i] <- delta*trapz(rfine[,i]^2)
}

C      <- 1 + intcovDhSy/mean(intysqr)
MS.amp <- C*mean(intrsqr)
MS.pha <- C*intetasqr - intmusqr
RSQR   <- MS.pha/(MS.amp+MS.pha)

return(list("MS.amp" = MS.amp, "MS.pha" = MS.pha, "RSQR" = RSQR, "C" = C))

}

trapz = function(x) {
n = length(x)
intx = sum(x) - 0.5*(x[1]+x[n])
return(intx)
}
```

Try the fda package in your browser

Any scripts or data that you put into this service are public.

fda documentation built on May 29, 2024, 11:26 a.m.