library(fda) library(ggplot2) # Acceleration for Handwrite Data # Set k-th sample k <- 5 rangval <- range(handwritTime) bbasis <- create.bspline.basis (rangval, norder = 6, nbasis = 33) BSmoothx <- smooth.basis(handwritTime, handwrit[, k, 1], bbasis) handwfdBx <- BSmoothx$fd fdBSmoothevalx <- eval.fd(fdobj = handwfdBx, evalarg = handwritTime) fdBSmoothevalder1x <- eval.fd(fdobj = deriv.fd(handwfdBx, 1), evalarg = handwritTime) fdBSmoothevalder2x <- eval.fd(fdobj = deriv.fd(handwfdBx, 2), evalarg = handwritTime) BSmoothy <- smooth.basis(handwritTime, handwrit[, k, 2], bbasis) handwfdBy <- BSmoothy$fd fdBSmoothevaly <- eval.fd(fdobj = handwfdBy, evalarg = handwritTime) fdBSmoothevalder1y <- eval.fd(fdobj = deriv.fd(handwfdBy, 1), evalarg = handwritTime) fdBSmoothevalder2y <- eval.fd(fdobj = deriv.fd(handwfdBy, 2), evalarg = handwritTime) # plot x11() plot(fdBSmoothevalx, fdBSmoothevaly, type = "l") # ggplot # dfB <- data.frame(time = handwritTime, sx = fdBSmoothevalx, sy = fdBSmoothevaly) # # x11() # p <- ggplot(dfB, aes(x = sx, y = sy)) # p <- p + geom_path(colour = 2) # p <- p + labs (x = "Position [mm]", y = "Position [mm]") # p <- p + theme_bw () # print (p) # Angular Acceleration # data angacc <- fdBSmoothevalder2x^2 + fdBSmoothevalder2y^2 # B-spline smoothing bbasisacc <- create.bspline.basis (rangval, norder = 4, nbasis = 60) BSmooth <- smooth.basis(handwritTime, angacc, bbasisacc) handwfdB <- BSmooth$fd fdBSmootheval <- eval.fd(fdobj = handwfdB, evalarg = handwritTime) # plot plotfit.fd(angacc, handwritTime, handwfdB, ylim = c(-1e-11,1e-10), cex.pch = 0.5) # ggplot # accdf <- data.frame(time = handwritTime, angacc = angacc) # # x11() # p <- ggplot(accdf, aes(x = time, y = angacc)) # p <- p + geom_point(colour = 2, size = 0.5) # p <- p + labs (x = "Time [ms]", y = "Acceleration [rad/ms^2]") # p <- p + geom_line(aes(y = fdBSmootheval), colour = 4) # p <- p + geom_hline(yintercept = 0) # p <- p + theme_bw () + scale_y_continuous(limits = c(-1e-11,1e-10)) # print (p) # positive smoothing ang <- angacc*10^10 bbasisacc <- create.bspline.basis (rangval, norder = 4, nbasis = 40) posBSmooth <- smooth.pos(handwritTime, ang, bbasisacc) poshandwfdB <- posBSmooth$Wfdobj posfdBSmootheval <- eval.posfd(Wfdobj = poshandwfdB, evalarg = handwritTime)/10^10 # plot # plotfit.fd doesn't work! plot(handwritTime, angacc, type = "p", pch = 20, ylim = c(-1e-11,1e-10)) lines(handwritTime, posfdBSmootheval, col = "red") # ggplot # x11() # p <- ggplot(accdf, aes(x = time, y = angacc)) # p <- p + geom_point(colour = 2, size = 0.5) # p <- p + labs (x = "Time [ms]", y = "Acceleration [rad/ms^2]") # p <- p + geom_line(aes(y = posfdBSmootheval), colour = 4) # p <- p + geom_hline(yintercept = 0) # p <- p + theme_bw () + scale_y_continuous(limits = c(-1e-11,1e-10)) # print (p) #-------------------------------------------------------------- # Monotone Smoothing # infantGrowth data (tibia - holenni kost) plot(infantGrowth[,2], xlab = "Days", ylab = "Length", type = "p") days <- c(infantGrowth[, 1]) tibia <- c(infantGrowth[, 2]) # B-spline smoothing # breaks <- days[-c(1, 40)] rangval <- range(days) norder <- 4 bbasis <- create.bspline.basis (rangval, norder = norder, breaks = days) curv.Lfd <- int2Lfd(2) curv.fdPar <- fdPar(bbasis, curv.Lfd, 1e-1) BSmooth <- smooth.basis(days, tibia, curv.fdPar) # plot plotfit.fd(tibia, days, BSmooth$fd) # The first derivative BSmoothder1 <- deriv.fd(BSmooth$fd, 1) plot(BSmoothder1) # # ggplot # df <- data.frame(x = days, y = tibia) # # xx <- seq(1, 40, length.out = 100) # fdBSmootheval <- eval.fd(fdobj = BSmooth$fd, evalarg = xx) # fdBSmoothevalder1 <- eval.fd(fdobj = deriv.fd(BSmooth$fd), evalarg = xx) # # dfsm <- data.frame(xx = xx, fdBSmootheval = fdBSmootheval, # fdBSmoothevalder1 = c(fdBSmoothevalder1))#, pos = c((fdBSmoothevalder1 > 0) + 1)) # # x11() # p0 <- ggplot (data = df, aes (x = x, y = y)) # p0 <- p0 + geom_point (colour = 2) # p0 <- p0 + geom_line(data = dfsm, aes (x = xx, y = fdBSmootheval), colour = 4) # p0 <- p0 + labs (x = "Day", y = "Length [mm]") # p0 <- p0 + theme_bw () + ggtitle("Baby's tibia") # print (p0) # # # derivative # x11() # p0 <- ggplot (data = dfsm, aes (x = xx, y = fdBSmoothevalder1)) # p0 <- p0 + geom_line(colour = 4) # p0 <- p0 + geom_rect(aes(xmin = 0, xmax = 40, ymin = -0.2, ymax = 0), alpha = 0.01, fill = "white") # p0 <- p0 + geom_hline(yintercept = 0, linetype = "dashed") # p0 <- p0 + labs (x = "Day", y = "D Length [mm/day]") # p0 <- p0 + theme_bw () + ggtitle("Baby's tibia") # print (p0) # monotone smoothing monocurv.fdPar <- fdPar(bbasis, curv.Lfd, 1e-4) monoBSmooth <- smooth.monotone (days, tibia, monocurv.fdPar) monotibfdB <- monoBSmooth$Wfdobj # plot # plotfit.fd doesn't work! plot(days, tibia, type = "p", pch = 20) xx <- seq(1, 40, length.out = 100) beta <- monoBSmooth$beta fdBSmootheval <- beta[1] + beta[2]*eval.monfd(Wfdobj = monotibfdB, evalarg = xx) lines(xx, fdBSmootheval, col = "red") # The first derivative fdBSmoothevalder1 <- beta[2]*eval.monfd(Wfdobj = monotibfdB, evalarg = xx, Lfdobj = 1) plot(xx, fdBSmoothevalder1, type = "l") lines(c(xx[1], xx[100]), c(0, 0), lty = "dashed") # ggplot # monodfsm <- data.frame(xx = xx, fdBSmootheval = fdBSmootheval, # fdBSmoothevalder1 = c(fdBSmoothevalder1)) # # x11() # p0 <- ggplot (data = df, aes (x = x, y = y)) # p0 <- p0 + geom_point (colour = 2) # p0 <- p0 + geom_line(data = monodfsm, aes (x = xx, y = fdBSmootheval), colour = 4) # p0 <- p0 + labs (x = "Day", y = "Length [mm]") # p0 <- p0 + theme_bw () + ggtitle("Baby's tibia") # print (p0) # # # # derivative # x11() # p0 <- ggplot (data = monodfsm, aes (x = xx, y = fdBSmoothevalder1)) # p0 <- p0 + geom_line(colour = 4) # #p0 <- p0 + geom_rect(aes(xmin = 0, xmax = 40, ymin = -0.2, ymax = 0), alpha = 0.01, fill = "white") # p0 <- p0 + geom_hline(yintercept = 0, linetype = "dashed") # p0 <- p0 + labs (x = "Day", y = "D Length [mm/day]") # p0 <- p0 + theme_bw () + ggtitle("Baby's tibia") # print (p0) #-------------------------------------------------------------- # Density estimation # Canadian Weather data(CanadianWeather) temp <- CanadianWeather$dailyAv[,,1] precip <- CanadianWeather$dailyAv[,,2] daytime <- (1:365)-0.5 df <- data.frame(daytime = daytime, temp = c(temp), precip = c(precip), place = rep(CanadianWeather$place, each = 365), region = rep(CanadianWeather$region, each = 365)) dfp <- subset(df, place == "St. Johns") # set up range for density rangval <- range(dfp$precip) # set up basis for W(x) breaks <- seq(rangval[1], rangval[2], length.out = 29) norder <- 6 bbasis <- create.bspline.basis (rangval, norder = norder, breaks = breaks) nbasis <- bbasis$nbasis curv.Lfd <- int2Lfd(3) # pozor! - tady jiny zpusob tvoreni fdPar objektu, predchozi zpusob nefunguje # implicitne tvori matici koeficientu nbasis x nbasis, chceme nbasis x 1 Wfd <- fd(matrix(0, nbasis, 1), bbasis) WfdParobj <- fdPar(Wfd, curv.Lfd, lambda = 1e-2) BSmooth <- density.fd(dfp$precip, WfdParobj) # plot density xx <- seq(rangval[1], rangval[2], length.out = 200) wval <- eval.fd(xx, BSmooth$Wfdobj) pval <- exp(wval)/BSmooth$C plot(xx, pval, type="l") # ggplot # densdf <- data.frame(x = dfp$precip , y = rep(0, length(dfp$precip))) # densdfsm <- data.frame(xx = xx, fdBSmootheval = c(pval)) # # x11() # p0 <- ggplot (data = densdf, aes (x = x, y = y)) # p0 <- p0 + geom_jitter (colour = 2, height = 0.002, size = 0.5) # p0 <- p0 + geom_line(data = densdfsm, aes (x = xx, y = fdBSmootheval), colour = 4) # p0 <- p0 + labs (x = "Precipitation", y = "Density") # p0 <- p0 + theme_bw () + ggtitle("St. Johns") # print (p0) # ------------------ # Problems # ------------------ ### Monotone Smoothing #load("absorb.RData") absorb <- read.csv("absorb.csv") x <- absorb$time y <- as.matrix(absorb[, 2:5]) ### Density estimation # load("turany.RData") df.turany.monthly <- read.csv("turany.csv") precip <- df.turany.monthly$precipitation