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))


dfp <- subset(df, place == "St. Johns")
rangeval <- range(dfp$daytime)

### Cross - validation

k <- 9 # number of basis
fbasis <- create.fourier.basis(rangeval, nbasis = k)
fdobjSmooth <- smooth.basis(dfp$daytime, dfp$precip, fbasis)
phi <- t(fdobjSmooth$y2cMap)
phiinv <- solve(t(phi)%*%phi)

# Smoothing matrix S
S <- phi%*%phiinv%*%t(phi)

# Estimate: yhat = S*y
yhat <- S%*%dfp$precip

sii <- diag(S)
OCV <- mean(((dfp$precip - yhat)/(1-sii))^2)

# or simply
OCV <- fdobjSmooth$gcv

### Confidence bands

k <- 9
fbasis <- create.fourier.basis(rangeval, nbasis = k)
fdobjSmooth <- smooth.basis(dfp$daytime, dfp$precip, fbasis)
fdobjSmootheval <- eval.fd(fdobj = fdobjSmooth$fd, evalarg = dfp$daytime)

# variance estimate
eps <- fdobjSmootheval - dfp$precip
var(eps)
# or simply
(sigma2 <- fdobjSmooth$SSE/(365 - k))

# Smoothing matrix S
phi <- t(fdobjSmooth$y2cMap)
phiinv <- solve(t(phi)%*%phi)
S <- phi%*%phiinv%*%t(phi)

Varyhat <- sigma2*S%*%t(S)
Lband <- fdobjSmootheval - 1.96*sqrt(diag(Varyhat))
Uband <- fdobjSmootheval + 1.96*sqrt(diag(Varyhat))

plot(dfp$daytime, dfp$precip, type = "p")
lines(dfp$daytime, fdobjSmootheval)
lines(dfp$daytime, Uband, lty = "dashed")
lines(dfp$daytime, Lband, lty = "dashed")

# ggplot
# x11()
# p0 <- ggplot (data = dfp, aes (x = daytime, y = precip, color = place))
# p0 <- p0 + geom_point ()
# p0 <- p0 + geom_line (aes(y = fdobjSmootheval), colour = 4)
# p0 <- p0 + geom_line (aes(y = Lband), colour = 4, linetype = "dashed")
# p0 <- p0 + geom_line (aes(y = Uband), colour = 4, linetype = "dashed")
# p0 <- p0 + labs (x = "Days", y = "Precipitation")
# p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns")
# print (p0)


# ------------------
# Problems
# ------------------

# 1. Melanoma Data
?melanoma

x <- melanoma$year
y <- melanoma$incidence

# residuals
melres <- lm(incidence ~ ., data = melanoma)$residuals

# 2. Canadian Weather Data
?CanadianWeather