Производные B сплайнового базиса в R - PullRequest
0 голосов
/ 18 февраля 2019

Я не совсем понимаю, как получить первую и вторую производную ab сплайновой базы для FDA в R. Я предпринял попытки, изменив «порядок» ряда и используя код, который работал для базисной системы Фурье.

Любая помощь будет очень цениться

gaitHip = gait[, , 1]
(gaittime <- as.numeric(dimnames(gaitHip)[[1]])*20)
gaitrange <- c(0,20)


# -----------  set up the harmonic acceleration operator  ----------

#harmaccelLfd <- vec2Lfd(c(0, (2*pi/20)^2, 0), rangeval=gaitrange)

Bsplineoperator <- int2Lfd(max(0, 2))

#  Set up basis for representing gaitHip data. 

Bsplinebasis = create.bspline.basis(gaitrange, nbasis=20, norder=4)


#  -------------------------------------------------------------------
#                 Choose level of smoothing using
#          the generalized cross-validation criterion
#  -------------------------------------------------------------------

#  set up range of smoothing parameters in log_10 units

gaitLoglam <- seq(-4,0,.25)
nglam   <- length(gaitLoglam)

gaitSmoothStats <- array(NA, dim=c(nglam, 3),
                         dimnames=list(gaitLoglam, c("log10.lambda", "df", "gcv") ) )
gaitSmoothStats[, 1] <- gaitLoglam

#  loop through smoothing parameters

for (ilam in 1:nglam) {
  gaitSmooth <- smooth.basisPar(gaittime, gaitHip, Bsplinebasis,
                                Lfdobj=Bsplineoperator, lambda=10^gaitLoglam[ilam])
  gaitSmoothStats[ilam, "df"]  <- gaitSmooth$df
  gaitSmoothStats[ilam, "gcv"] <- sum(gaitSmooth$gcv)
  # note: gcv is a matrix in this case
}

#  display and plot GCV criterion and degrees of freedom

op <- par(mfrow=c(1,1))
plot(gaitLoglam, gaitSmoothStats[, "gcv"],pch=19, type="b",
     xlab=expression("Log"[10]*lambda), ylab="GCV Criterion",
     main="Gait Smoothing", log="y")

gaitfd <- smooth.basisPar(gaittime, gaitHip,
                          Bsplinebasis,Lfdobj=Bsplineoperator, lambda=1e-1)$fd
names(gaitfd$fdnames) <- c("Normalized time", "Child", "Angle")
gaitfd$fdnames[[3]] <- c("Hip")

str(gaitfd)

Derivative1=deriv(gaitfd, Lfdobj=int2Lfd(1))
Derivative2=deriv(gaitfd, Lfdobj=int2Lfd(2))


CoefficientsHip=gaitfd[["coefs"]]
CoefficientsHip1=Derivative1[["coefs"]]
CoefficientsHip2=Derivative2[["coefs"]]

gaitHip <- структура (c (37, 36, 33, 29, 23, 18, 15, 12, 9, 6, 4, 6, 13, 22, 31, 38, 43, 44, 40, 35, 47, 46, 42, 34, 27, 21, 16, 12, 4, -1, -2, 2, 12, 24, 34, 42, 47, 48, 45, 43, 46, 44, 39, 34, 33, 27, 21, 15, 7, 1, -3, -1, 10, 22, 34, 43, 51, 53, 50, 49, 37, 36, 27, 20, 15, 15, 10, 6, 3, 0, 3, 4, 10, 14, 23, 31, 33, 35, 34, 32, 20, 18, 11, 8, 7, 5, 1, 0, -2, -4, -4, 1, 10, 18, 27, 36, 34, 33, 32, 28, 57, 48, 44, 35, 31, 27, 21,18, 14, 10, 8, 8, 11, 21, 32, 40, 49, 55, 56, 55, 46, 38, 33, 25, 18, 15, 14, 10, 8, 4, 0, 3,9, 19, 31, 37, 46, 45, 46, 41, 46, 46, 43, 40, 36, 30, 24, 17, 12, 7, 1, 3, 12, 24, 35, 44, 49,49, 47, 43, 46, 42, 37, 34, 31, 25, 19, 17, 13, 10, 5, 7, 13, 21, 29, 39, 47, 51, 52, 47, 35, 34,29, 28, 19, 15, 13, 6, 0, -1, 2, 11, 20, 30, 39, 45, 48, 48, 43, 39, 38, 37, 33, 29, 26, 20, 14, 9, 3, -2, -1, 5, 15, 24, 36, 44, 47, 45, 42, 40, 35, 31, 29, 26, 22, 19, 15, 11, 6, 3, 0, 4, 16, 27, 35, 41, 44, 45, 40, 39, 34, 31, 27, 23, 19, 15, 13, 10, 6, 1, 0, 3, 8, 16, 27, 34, 40, 41, 39, 36, 43, 41, 36, 31, 26, 20, 16, 13, 10, 8, 4, 12, 19, 27, 37, 41, 44, 44, 41, 37, 43, 37, 35, 28, 26, 21, 20, 15, 13, 7, 3, 8, 16, 25, 35, 43, 49, 50, 45, 46, 40, 41, 36, 32, 27, 20, 17, 10, 5, -2, -4, -2, 7, 24, 38, 49, 57, 59, 54, 46, 51, 49, 45, 39, 31, 23, 16, 9, 4, 0, -5, 1, 10, 21, 33, 44, 51, 55, 56, 51, 52, 46, 41, 35, 31, 24, 18, 12, 6, 3, 2, 5, 12, 22, 32, 42, 46, 48, 49, 46, 36, 33, 28, 22, 18, 13, 11, 7, 1, -2, -3, -3, 2, 13, 22, 32, 39, 41, 38, 34, 35, 37, 33,27, 22, 14, 9, 6, 4, -1, 1, 2, 11, 17, 34, 35, 40, 43, 43, 42, 46, 38, 30, 23, 17, 13, 7, 3, 0, -3, -3, 1, 11, 20, 36, 48, 55, 57, 56, 50, 43, 41, 37, 30, 24, 16, 12, 8, 5, 2, 1, 8, 16, 24, 33, 42, 48, 48, 48, 46, 55, 51, 47, 41, 35, 30, 26, 23, 20, 11, 8, 10, 19, 35, 41, 52, 57, 61, 63, 58, 39, 38, 31, 27, 21, 14, 9, 9, 6, 4, 2, 3, 13, 30, 37, 47, 53, 53, 49, 44, 37, 34, 30, 27, 26, 19, 15, 10, 4, 0, 0, 4, 12, 22, 31, 40, 43, 43, 38, 38, 36, 33, 28, 22, 18, 13, 11, 7, 1, -2, -3, -3, 2, 13, 22, 32, 39, 41, 38, 34, 36, 33,30, 28, 21, 15, 8, 1, -5, -11, -12, -7, 4, 16, 26, 37, 44, 47, 44, 37, 42, 40, 40, 34, 23,15, 11, 7, 5, 6, 8, 12, 22, 33, 43, 51, 57, 58, 54, 46, 38, 34, 30, 23, 17, 12, 8, 4, 1, 0,-4, -4, 1, 10, 22, 32, 38, 41, 41, 40, 46, 47, 44, 37, 29, 23, 19, 14, 8, 3, -2, 1, 12, 26, 39, 48, 52, 48, 43, 42, 54, 48, 44, 37, 30, 27, 21, 18, 15, 11, 12, 20, 30, 41, 50, 56, 61, 59, 57, 58, 52, 44, 44, 33, 28, 27, 23, 24, 19, 15, 15, 16, 24, 32, 43, 52, 58, 59, 57, 52, 32, 28, 26, 22, 19, 13, 8, 5, -1, -6, -5, 0, 12, 22, 30, 36, 39, 36, 30, 29, 46, 41, 38, 31, 25, 20, 13,7, 1, -4, -3, -2, 8, 22, 34, 45, 53, 57, 55, 43, 46, 44, 40, 35, 31, 25, 19, 15, 10, 5, 3, 6, 14, 27, 40, 48, 53, 53, 50, 47, 48, 42, 42, 35, 30, 23, 19, 14, 9, 4, -1, 3, 9, 25, 37,45, 52, 53, 52, 46, 44, 41, 38, 32, 24, 18, 10, 6, 3, 0, -2, -1, 6, 20, 36, 44, 49, 46, 38,35, 55, 56, 51, 46, 41, 36, 30, 25, 21, 15, 8, 5, 9, 19, 31, 43, 52, 56, 59, 59, 48, 50, 47, 42, 37, 29, 22, 14, 8, 5, 8, 15, 24, 36, 51, 59, 63, 64, 61, 55), .Dim = c(20L, 39L), .Dimnames = list (c («0,025», «0,075», «0,125», «0,175», «0,225», «0,275», «0,325», «0,375», «0,425», «0,475 "," 0,525 "," 0,575 "," 0,625 "," 0,675 "," 0,725 "," 0,775 "," 0,825 "," 0,875 "," 0,925 "," 0,975 "), c (" boy1 ",«boy2», «boy3», «boy4», «boy5», «boy6», «boy7», «boy8», «boy9», «boy10», «boy11», «boy12», «boy13», «boy14»"," boy15 "," boy16 "," boy17 "," boy18 "," boy19 "," boy20 "," boy21 "," boy22 "," boy23 "," boy24 "," boy25 "," boy26 ",«boy27», «boy28», «boy29», «boy30», «boy31», «boy32», «boy33», «boy34», «boy35», «boy36», «boy37», «boy38», «boy39»«))) </p>

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...