library(fda) library(ggplot2) # 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)) # Select just St. Johns dfp <- subset(df, place == "St. Johns") ### Plot these data x11() plot(dfp$daytime, dfp$temp) # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # ------------------ # The monomial basis # ------------------ # for n = 1 rangeval <- range(dfp$daytime) mbasis <- create.monomial.basis(rangeval, nbasis = 1) xx <- dfp$daytime[1:100] # Plot basis functions plot(mbasis) monoSmooth <- smooth.basis(dfp$daytime, dfp$temp, mbasis) # plot data and smoothed function plot(dfp$daytime, dfp$temp) lines(monoSmooth$fd) # fdmonoSmootheval <- eval.fd(fdobj = monoSmooth$fd, evalarg = dfp$daytime) # fdmono <- eval.basis(basisobj = mbasis, evalarg = xx) # # ebchan <- fdmono*matrix(rep(monoSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdmono), each = length(xx))) # # # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Monomial Basis") # print (p0) # # # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdmonoSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # # for n = 2 mbasis <- create.monomial.basis(rangeval, nbasis = 2) plot(mbasis) monoSmooth <- smooth.basis(dfp$daytime, dfp$temp, mbasis) plot(dfp$daytime, dfp$temp) lines(monoSmooth$fd) # fdmonoSmootheval <- eval.fd(fdobj = monoSmooth$fd, evalarg = dfp$daytime) # fdmono <- eval.basis(basisobj = mbasis, evalarg = xx) # ebchan <- fdmono*matrix(rep(monoSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdmono), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Monomial Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdmonoSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # for n = 3 mbasis <- create.monomial.basis(rangeval, nbasis = 3) plot(mbasis) monoSmooth <- smooth.basis(dfp$daytime, dfp$temp, mbasis) plot(dfp$daytime, dfp$temp) lines(monoSmooth$fd) # fdmonoSmootheval <- eval.fd(fdobj = monoSmooth$fd, evalarg = dfp$daytime) # fdmono <- eval.basis(basisobj = mbasis, evalarg = xx) # ebchan <- fdmono*matrix(rep(monoSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdmono), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Monomial Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdmonoSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # for n = 6 mbasis <- create.monomial.basis(rangeval, nbasis = 6) plot(mbasis) monoSmooth <- smooth.basis(dfp$daytime, dfp$temp, mbasis) plot(dfp$daytime, dfp$temp) lines(monoSmooth$fd) # fdmonoSmootheval <- eval.fd(fdobj = monoSmooth$fd, evalarg = dfp$daytime) # fdmono <- eval.basis(basisobj = mbasis, evalarg = xx) # ebchan <- fdmono*matrix(rep(monoSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdmono), each = length(xx))) # # # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Monomial Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdmonoSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # First derivative for mono 6 mono1d <- deriv.fd(monoSmooth$fd, 1) plot(mono1d) # fdmonoSmootheval <- eval.fd(fdobj = monoSmooth$fd, evalarg = dfp$daytime) # fdobjDer1 <- eval.fd(fdobj = deriv.fd(monoSmooth$fd, 1), evalarg = dfp$daytime) # fdobjDer2 <- eval.fd(fdobj = deriv.fd(monoSmooth$fd, 2), evalarg = dfp$daytime) # # dfp$tempder1mono <- fdobjDer1 # dfp$tempder2mono <- fdobjDer2 # # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = tempder1mono)) # p0 <- p0 + geom_line (colour = 4) # p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed") # p0 <- p0 + labs (x = "Days", y = "D Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # Second derivative for mono 6 mono2d <- deriv.fd(monoSmooth$fd, 2) plot(mono2d) # p0 <- ggplot (data = dfp, aes (x = daytime, y = tempder2mono)) # p0 <- p0 + geom_line (colour = 4) # p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed") # p0 <- p0 + labs (x = "Days", y = "D2 Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # ----------------- # The Fourier basis # ----------------- # for n = 1 fbasis <- create.fourier.basis(rangeval, nbasis = 1, period = 365) xx <- dfp$daytime plot(fbasis) FSmooth <- smooth.basis(dfp$daytime, dfp$temp, fbasis) plot(dfp$daytime, dfp$temp) lines(FSmooth$fd) # FSmooth$fd <- FSmooth$fd # fdFSmootheval <- eval.fd(fdobj = FSmooth$fd, evalarg = dfp$daytime) # fdF <- eval.basis(basisobj = fbasis, evalarg = xx) # # ebchan <- fdF*matrix(rep(FSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdF), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Fourier Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdFSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # for n = 2 fbasis <- create.fourier.basis(rangeval, nbasis = 3, period = 365) plot(fbasis) FSmooth <- smooth.basis(dfp$daytime, dfp$temp, fbasis) plot(dfp$daytime, dfp$temp) lines(FSmooth$fd) # fdFSmootheval <- eval.fd(fdobj = FSmooth$fd, evalarg = dfp$daytime) # fdF <- eval.basis(basisobj = fbasis, evalarg = xx) # # ebchan <- fdF*matrix(rep(FSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdF), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Fourier Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdFSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # for n = 3 fbasis <- create.fourier.basis(rangeval, nbasis = 5, period = 365) plot(fbasis) FSmooth <- smooth.basis(dfp$daytime, dfp$temp, fbasis) plot(dfp$daytime, dfp$temp) lines(FSmooth$fd) # fdFSmootheval <- eval.fd(fdobj = FSmooth$fd, evalarg = dfp$daytime) # fdF <- eval.basis(basisobj = fbasis, evalarg = xx) # # ebchan <- fdF*matrix(rep(FSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdF), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Fourier Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdFSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # for n = 8 fbasis <- create.fourier.basis(rangeval, nbasis = 15, period = 365) plot(fbasis) FSmooth <- smooth.basis(dfp$daytime, dfp$temp, fbasis) plot(dfp$daytime, dfp$temp) lines(FSmooth$fd) # fdFSmootheval <- eval.fd(fdobj = FSmooth$fd, evalarg = dfp$daytime) # fdF <- eval.basis(basisobj = fbasis, evalarg = xx) # # ebchan <- fdF*matrix(rep(FSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdF), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Fourier Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdFSmootheval), colour = 4) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # First derivative for Four 8 F1d <- deriv.fd(FSmooth$fd, 1) plot(F1d) # fdFSmootheval <- eval.fd(fdobj = FSmooth$fd, evalarg = dfp$daytime) # fdobjDer1 <- eval.fd(fdobj = deriv.fd(FSmooth$fd, 1), evalarg = dfp$daytime) # fdobjDer2 <- eval.fd(fdobj = deriv.fd(FSmooth$fd, 2), evalarg = dfp$daytime) # # dfp$tempder1F <- fdobjDer1 # dfp$tempder2F <- fdobjDer2 # # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = tempder1F)) # p0 <- p0 + geom_line (colour = 2, linetype = "dashed") # p0 <- p0 + geom_line (aes(y = tempder1), colour = 4) # p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed") # p0 <- p0 + labs (x = "Days", y = "D Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # First derivative for Four 8 F2d <- deriv.fd(FSmooth$fd, 2) plot(F2d) # p0 <- ggplot (data = dfp, aes (x = daytime, y = tempder2F)) # p0 <- p0 + geom_line (colour = 2, linetype = "dashed") # p0 <- p0 + geom_line (aes(y = tempder2), colour = 4) # p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed") # p0 <- p0 + labs (x = "Days", y = "D2 Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # ----------------- # The Bspline basis # ----------------- # As the first, 3 knots and local constants breaks <- cumsum(c(0, 212, 153)) xx <- dfp$daytime # for n = 1 norder <- 1 # This means I have 1+1 basis functions nbasis <- 2 bbasis <- create.bspline.basis (rangeval, norder = norder, breaks = breaks) plot(bbasis) BSmooth <- smooth.basis(dfp$daytime, dfp$temp, bbasis) plot(dfp$daytime, dfp$temp) lines(BSmooth$fd) # fdBSmootheval <- eval.fd(fdobj = BSmooth$fd, evalarg = dfp$daytime) # fdB <- eval.basis(basisobj = bbasis, evalarg = xx) # # ebchan <- fdB*matrix(rep(BSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdB), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Bsplines Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdBSmootheval), colour = 4) # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # As the second, knots monthly and local constants breaks <- cumsum(c(0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)) xx <- dfp$daytime # m = 1 norder <- 1 # This means I have 11+1 basis functions nbasis <- 12 bbasis <- create.bspline.basis (rangeval, norder = norder, breaks = breaks) plot(bbasis) BSmooth <- smooth.basis(dfp$daytime, dfp$temp, bbasis) plot(dfp$daytime, dfp$temp) lines(BSmooth$fd) # BSmooth$fd <- BSmooth$fd # fdBSmootheval <- eval.fd(fdobj = BSmooth$fd, evalarg = dfp$daytime) # fdB <- eval.basis(basisobj = bbasis, evalarg = xx) # # ebchan <- fdB*matrix(rep(BSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdB), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Bsplines Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdBSmootheval), colour = 4) # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # m = 2 norder <- 2 # This means I have 11+2 basis functions nbasis <- 13 bbasis <- create.bspline.basis (rangeval, norder = norder, breaks = breaks) plot(bbasis) BSmooth <- smooth.basis(dfp$daytime, dfp$temp, bbasis) plot(dfp$daytime, dfp$temp) lines(BSmooth$fd) # fdBSmootheval <- eval.fd(fdobj = BSmooth$fd, evalarg = dfp$daytime) # fdB <- eval.basis(basisobj = bbasis, evalarg = xx) # # ebchan <- fdB*matrix(rep(BSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdB), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Bsplines Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdBSmootheval), colour = 4) # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # m = 3 norder <- 3 # This means I have 11+3 basis functions nbasis <- 14 bbasis <- create.bspline.basis (rangeval, norder = norder, breaks = breaks) plot(bbasis) BSmooth <- smooth.basis(dfp$daytime, dfp$temp, bbasis) plot(dfp$daytime, dfp$temp) lines(BSmooth$fd) # fdBSmootheval <- eval.fd(fdobj = BSmooth$fd, evalarg = dfp$daytime) # fdB <- eval.basis(basisobj = bbasis, evalarg = xx) # # ebchan <- fdB*matrix(rep(BSmooth$fd$coefs, each = length(xx)), nrow = length(xx)) # basisdf <- data.frame(bs = c(ebchan), x = xx, basis = rep(colnames(fdB), each = length(xx))) # samotna baze # x11() # p0 <- ggplot (data = basisdf, aes (x = x, y = bs, colour = basis)) # p0 <- p0 + geom_line () # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Basis") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("Bsplines Basis") # print (p0) # vyhlazena data # x11() # p0 <- ggplot (data = dfp, aes (x = daytime, y = temp, color = place)) # p0 <- p0 + geom_point () # p0 <- p0 + geom_line (aes(y = fdBSmootheval), colour = 4) # p0 <- p0 + geom_vline (xintercept = breaks, linetype = "dotted", size = 0.3) # p0 <- p0 + labs (x = "Days", y = "Temperature") # p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns") # print (p0) # ------------------ # Problems # ------------------ # 1. Canadian Weather Data ?CanadianWeather temp <- CanadianWeather$dailyAv[,,1] daytime <- (1:365)-0.5 # set knots breaks <- c(cumsum(c(0, 82, 91, 91, 91)), 365) # set knots breaks <- seq(0,365,5) #2. Refinery Data ?refinery # if doesn't work: source("refinery.R") x <- refinery$Time y <- refinery$Tray47 # set knots breaks <- c(0, 33, 67, 98, 130, 162, 193) # set knots - double 67 breaks <- c(0, 33, 67, 67, 98, 130, 162, 193) # set knots - triple 67 breaks <- c(0, 33, 67, 67, 67, 98, 130, 162, 193)