Skip to content
Snippets Groups Projects
Commit 103aa8ec authored by pbac's avatar pbac
Browse files

Different updates, on the way to v0.9.1

parent 5a9f9ad5
Branches
No related tags found
No related merge requests found
......@@ -14,6 +14,7 @@ Imports:
Rcpp (>= 0.12.18),
R6 (>= 2.2.2),
splines (>= 3.1.1),
pbs,
digest,
LinkingTo: Rcpp, RcppArmadillo
Suggests:
......
......@@ -48,7 +48,7 @@ as.data.list <- function(object){
#' as.data.list(X)
#'
#' # Convert a dataframe with time, forecast and an observed variable
#' X <- data.frame(t=1:10, x.k1=1:10, x.k2=10:1, y.obs=1:10, y.k1=1:10, y.k2=1:10)
#' X <- data.frame(t=1:10, x.k1=1:10, x.k2=10:1, yobs=1:10, y.k1=1:10, y.k2=1:10)
#' as.data.list(X)
#'
#' # Can be converted back and forth
......
......@@ -17,18 +17,19 @@
#'
#' @param X data.frame (as part of data.list) with horizons as columns named \code{kxx} (i.e. one for each horizon)
#' @param Boundary.knots The value is NA: then the boundaries are set to the range of each horizons (columns in X). See \code{?splines::bs}
#' @param intercept Default value is TRUE: in an onlineforecast model there is no intercept per default (set by \code{one()}). It's default to FALSE in \code{?splines::bs}.
#' @param intercept See \code{?splines::bs}.
#' @param df See \code{?splines::bs}
#' @param knots See \code{?splines::bs}
#' @param degree See \code{?splines::bs}
#' @param bknots Is just a short for Boundary.knots and replace Boundary.knots (if Boundary.knots is not given)
#' @param periodic Default FALSE. If TRUE, then \code{pbs::pbs} is called and periodic splines are generated.
#' @return List of data frames with the computed base splines, each with columns for the same horizons as in X
#' @rdname bs
#' @examples
#'
#' # How to make a diurnal curve using splines
#' # Select first 54 hours from the load data
#' D <- subset(Dbuilding, 1:54, kseq=1:4)
#' D <- subset(Dbuilding, 1:76, kseq=1:4)
#' # Make the hour of the day as a forecast input
#' D$tday <- make_tday(D$t, kseq=1:4)
#' D$tday
......@@ -52,9 +53,20 @@
#' # Such that at the transform stage will give the same as above
#' model$transform_data(D)
#'
#' # Periodic splines are useful for modelling a diurnal harmonical functions
#' L <- bspline(D$tday, bknots=c(0,24), df=4, periodic=TRUE)
#' # or
#' L <- pbspline(D$tday, bknots=c(0,24), df=4)
#' # Note, how it has to have high enough df, else it generates an error
#'
#' # Plot
#' plot(D$t, L$bs1$k1, type="s")
#' for(i in 2:length(L)){
#' lines(D$t, L[[i]]$k1, col=i, type="s")
#' }
#'
#' @export
bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots = NULL, degree = 3, bknots = NA) {
bspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots = NULL, degree = 3, bknots = NA, periodic = FALSE) {
# bknots is just a short for Boundary.knots and replace if Boundary.knots are not given.
if(is.na(Boundary.knots)){
Boundary.knots <- bknots
......@@ -64,7 +76,7 @@ bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots =
# Call again for each element
val <- lapply(1:length(X), function(i) {
bspline(X[[i]], df = df, knots = knots, degree = degree, intercept = intercept,
Boundary.knots = Boundary.knots)
Boundary.knots = Boundary.knots, periodic = periodic)
})
nams(val) <- nams(X)
return(val)
......@@ -77,8 +89,21 @@ bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots =
if (is.na(Boundary.knots[1])) {
Boundary.knots <- range(X[, i], na.rm=TRUE)
}
if(periodic){
# Periodic splines (annoyingly it use 'missing()' to conclude if the argument was given, not 'is.null()' as bs does)
if(is.null(knots)){
# Call without knots
spls <- pbs::pbs(X[, i], Boundary.knots = Boundary.knots, degree = degree, df = df,
intercept = intercept)
}else{
# Call without df
spls <- pbs::pbs(X[, i], Boundary.knots = Boundary.knots, degree = degree,
knots = knots, intercept = intercept)
}
}else{
spls <- splines::bs(X[, i], Boundary.knots = Boundary.knots, degree = degree, df = df,
knots = knots, intercept = intercept)
}
return(spls)
})
# Now we have a bs value in L for each horizon
......@@ -95,3 +120,24 @@ bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots =
nams(L) <- pst("bs", 1:length(L))
return(L)
}
#' Wrapper for \code{bspline} with \code{periodic=TRUE}
#'
#' Simply a wrapper.
#'
#' @title Wrapper for \code{bspline} with \code{periodic=TRUE}
#' @param X see \code{?bspline}
#' @param Boundary.knots see \code{?bspline}
#' @param intercept see \code{?bspline}
#' @param df see \code{?bspline}
#' @param knots see \code{?bspline}
#' @param degree see \code{?bspline}
#' @param bknots see \code{?bspline}
#' @family Transform stage functions
#'
#' @export
pbspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots = NULL, degree = 3, bknots = NA){
bspline(X, Boundary.knots = Boundary.knots, degree = degree, df = df,
knots = knots, intercept = intercept, bknots = bknots, periodic = TRUE)
}
......@@ -11,8 +11,8 @@
#' \item{t}{Time in GMT as POSIXct}
#' \item{heatload}{The heatload of a single family building in W}
#' \item{heatloadtotal}{The average heatload of a 16 single family buildings in W}
#' \item{Ta.obs}{Observed ambient temperature at the weather station in Celcius}
#' \item{I.obs}{Observed global radiation at the weather station in W/m^2}
#' \item{Taobs}{Observed ambient temperature at the weather station in Celcius}
#' \item{Iobs}{Observed global radiation at the weather station in W/m^2}
#' \item{Ta}{Weather forecasts of ambient temperature up to 36 hours ahead from DMI in Celcius}
#' \item{Ta}{Weather forecasts of global radiation up to 36 hours ahead from DMI in W/m^2}
#' }
......
......@@ -25,12 +25,12 @@
#' # The time vector
#' time <- seq(ct("2019-01-01"),ct("2019-01-02"),by=3600)
#' # Observations time series (as vector)
#' x.obs <- rnorm(length(time))
#' xobs <- rnorm(length(time))
#' # Forecast input as data.frame
#' X <- data.frame(matrix(rnorm(length(time)*3), ncol=3))
#' names(X) <- pst("k",1:3)
#'
#' D <- data.list(t=time, x.obs=x.obs, X=X)
#' D <- data.list(t=time, xobs=xobs, X=X)
#'
#' # Check it
#' check(D)
......@@ -58,7 +58,7 @@ data.list <- function(...) {
#' # Use the data.list with building heat load
#' D <- Dbuilding
#' # Take a subset for the example
#' D <- subset(D, 1:10, nms=c("t","Ta.obs","Ta","I.obs","I"), kseq=1:3)
#' D <- subset(D, 1:10, nms=c("t","Taobs","Ta","Iobs","I"), kseq=1:3)
#'
#' # Take subset index 2:4
#' subset(D, 2:4)
......@@ -74,7 +74,7 @@ data.list <- function(...) {
#' subset(D, nms=c("I","Ta"), kseq = 1)
#'
#' # Lag the forecasts such that they are aligned in time with observations
#' subset(D, nms=c("Ta.obs","Ta"), kseq = 2:3, lagforecasts = TRUE)
#' subset(D, nms=c("Taobs","Ta"), kseq = 2:3, lagforecasts = TRUE)
#'
#' # The order follows the order in nms
#' subset(D, nms=c("Ta","I"), kseq = 2)
......@@ -87,10 +87,10 @@ data.list <- function(...) {
#'
#' # A scatter plot between the forecast and the observations
#' # (try lagforecasts = FALSE and see the difference)
#' plot(X$Ta$k10, X$Ta.obs)
#' plot(X$Ta$k10, X$Taobs)
#'
#' # Fit a model for the 10-step horizon
#' abline(lm(Ta.obs ~ Ta.k10, X), col=2)
#' abline(lm(Taobs ~ Ta.k10, X), col=2)
#'
#' @export
subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts = FALSE, pattern = NA, ...) {
......@@ -200,7 +200,7 @@ subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
#' #' # Use the data.list with building heat load
#' D <- Dbuilding
#' # Take a subset
#' D <- subset(D, 1:5, nms=c("t","Ta.obs","Ta","I.obs","I"), kseq=1:3)
#' D <- subset(D, 1:5, nms=c("t","Taobs","Ta","Iobs","I"), kseq=1:3)
#'
#' # Convert to a data.frame, note the names of the forecasts are appended .kxx (i.e. for Ta and I)
#' as.data.frame(D)
......
......@@ -44,7 +44,7 @@
#' # Make two plots (and set the space for the legend)
#' plot_ts(D, c("heatload","Ta"), kseq=c(1,24), legendspace=11)
#' # Only the Ta observations
#' plot_ts(D, c("heatload","Ta.obs$"), kseq=c(1,24), legendspace=11)
#' plot_ts(D, c("heatload","Taobs$"), kseq=c(1,24), legendspace=11)
#'
#' # Give labels
#' plot_ts(D, c("heatload","Ta"), kseq=c(1,24), xlab="Time", ylabs=c("Heat (kW)","Temperature (C)"))
......
......@@ -20,7 +20,7 @@
#' D <- Dbuilding
#' plotly_ts(D, c("heatload","Ta"), kseq=c(1,24))
#' plotly_ts(D, c("heatload","Ta"), kseq=c(1,24))
#' plotly_ts(D, c("heatload","Ta$|Ta.obs$"), kseq=c(1,24))
#' plotly_ts(D, c("heatload","Ta$|Taobs$"), kseq=c(1,24))
#' }
#'
#' @export
......
......@@ -202,7 +202,7 @@ rls_fit <- function(prm=NA, model, data, scorefun = NA, returnanalysis = TRUE,
# The estimated coefficients
Lfitval <- getse(Lresult, "Theta", fun=as.data.frame)
# Return the model validation data
invisible(structure(list(Yhat = Yhat, model = model$clone_deep(), data = data, Lfitval = Lfitval, scoreval = scoreval), class = c("forecastmodel_fit","rls_fit")))
invisible(structure(list(Yhat = Yhat, model = model$clone_deep(), data = data, datatr = datatr, Lfitval = Lfitval, scoreval = scoreval), class = c("forecastmodel_fit","rls_fit")))
}else{
# Only the summed score returned
val <- sum(scoreval, na.rm = TRUE)
......
......@@ -133,7 +133,7 @@ rls_update <- function(model, datatr = NA, y = NA, runcpp=TRUE) {
}
#
# Give names to the matrices
colnames(P) <- names(datatr)
colnames(R) <- names(datatr)
colnames(Theta) <- names(datatr)
# Return the fit and result
return(list(fit=list(k=k, theta=theta, R=R, yhat=yhat[(n-kmax+1):n]),
......
......@@ -13,7 +13,7 @@
#' @title Computes the RMSE score.
#' @param x a numerical vector of residuals.
#' @return The RMSE score.
#' @seealso \code{\link{score_for_k}()} for calculation of a score for the k'th horizon, and \code{\link{score_fit}()} which takes a forecastmodel fit and returns score taking scoreperiod etc. into account.
#' @seealso \code{\link{score}()} for calculation of a score for the k'th horizon, and \code{\link{score_fit}()} which takes a forecastmodel fit and returns score taking scoreperiod etc. into account.
#' @name rmse
#' @examples
#'
......
......@@ -8,7 +8,7 @@
#' @param scoreperiod as an index (logical or integer) defining which points to inlude in the score calculation. If NA, the \code{scoreperiod} from the \code{fit$model} object is used.
#' @param usecomplete Only use points where
#' @param scorefun The score function applied, per default \code{\link{rmse}}.
#' @seealso \code{\link{rmse}} and \code{\link{score_for_k}} which are used in this function.
#' @seealso \code{\link{rmse}} and \code{\link{score}} which are used in this function.
#' @return A list with:
#' - \code{scoreval} is the score value
#' - \code{scoreperiod} is the score period used (can be different that the one in arguments \code{fit} or \code{scoreperiod})
......@@ -41,7 +41,7 @@ score_fit <- function(fit, scoreperiod = NA, usecomplete = TRUE, scorefun = rmse
}
# Calculate the score
tmp <- score_for_k(Residuals, scoreperiod, usecomplete)
tmp <- score(Residuals, scoreperiod, usecomplete)
scoreval <- tmp$scoreval
scoreperiod <- tmp$scoreperiod
......
......@@ -2,7 +2,7 @@
#library(devtools)
#document()
#load_all(as.package("../../onlineforecast"))
#?score_for_k
#?score
#' Calculates the score for each horizon for a matrix with residuals for each horizon.
#'
......@@ -24,15 +24,15 @@
#' Resid <- residuals(Yhat, y)
#'
#' # Calculate the score for the k1 horizon
#' score_for_k(Resid)$scoreval
#' score(Resid)$val
#'
#' # The first values were excluded, since there are NAs
#' head(Resid)
#' score_for_k(Resid)$scoreperiod
#' score(Resid)$scoreperiod
#'
#' @importFrom stats complete.cases
#' @export
score_for_k <- function(Residuals, scoreperiod = NA, usecomplete = TRUE, scorefun = rmse){
score <- function(Residuals, scoreperiod = NA, usecomplete = TRUE, scorefun = rmse){
# If no scoreperiod is given, then use all
if(is.na(scoreperiod[1])){
scoreperiod <- rep(TRUE,nrow(Residuals))
......@@ -50,10 +50,10 @@ score_for_k <- function(Residuals, scoreperiod = NA, usecomplete = TRUE, scorefu
scoreperiod <- scoreperiod & complete.cases(Residuals)
}
# Calculate the objective function for each horizon
scoreval <- sapply(1:ncol(Residuals), function(i){
val <- sapply(1:ncol(Residuals), function(i){
scorefun(Residuals[scoreperiod,i])
})
nams(scoreval) <- gsub("h","k",nams(Residuals))
nams(val) <- gsub("h","k",nams(Residuals))
#
return(list(scoreval=scoreval,scoreperiod=scoreperiod))
return(list(val=val,scoreperiod=scoreperiod))
}
No preview for this file type
......@@ -53,6 +53,8 @@ names(data[["Heatload"]]) <- pst("house", 1:16)
data[["cosAoi"]] <- data_or[, "cosAoi.obs"]
data[["sunElevation"]] <- data_or[, "sunElevation.obs"]
# .obs to obs
names(data) <- gsub("\\.","",names(data))
# # The time of day
# ncol <- ncol(data$Ta)
......@@ -78,9 +80,9 @@ data$totaln <- sapply(1:nrow(data$Heatload), function(i){
plot(data$t, data$totaln)
# Write for building heat load forecasting
#Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","Heatload","Ta","I","Ws","Wd","Ta.obs","I.obs","Wd.obs","Ws.obs","cosAoi","sunElevation","tday"))
#Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","Heatload","Ta","I","Ws","Wd","Taobs","Iobs","Wdobs","Wsobs","cosAoi","sunElevation","tday"))
data$heatload <- data$Heatload$house9
Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","heatload","heatloadtotal","Ta.obs","I.obs","Ta","I"))
Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","heatload","heatloadtotal","Taobs","Iobs","Ta","I"))
rownames(Dbuilding$Ta) <- NULL
Dbuilding$Ta <- Dbuilding$Ta[ ,1:36]
rownames(Dbuilding$I) <- NULL
......
......@@ -16,7 +16,6 @@
## install.packages("plotly")
#----------------------------------------------------------------
# Go
library(devtools)
......@@ -67,6 +66,14 @@ build(".", vignettes=TRUE)
install.packages("../onlineforecast_0.9.1.tar.gz")
library(onlineforecast)
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# Build binary package
#system("R CMD INSTALL --build ../onlineforecast_0.9.1.tar.gz")
# ----------------------------------------------------------------
# ----------------------------------------------------------------
......
......@@ -326,14 +326,14 @@ Now the residuals can be calculated and the score:
# Use the residuals function
R <- residuals(D$Yhat1, D$y)
# And the score as a function of the horizon
score_for_k(R, scoreperiod=ok)$scoreval
score(R, scoreperiod=ok)$val
```
Calculated the score (default is RMSE) for all models:
```{r}
RMSE <- sapply(nms, function(nm){
score_for_k(residuals(D[[nm]],D$y), ok)$scoreval
score(residuals(D[[nm]],D$y), ok)$val
})
```
......@@ -386,7 +386,7 @@ fittmp <- rls_fit(model$prm, model, D)
Finally, the score can be calculated on the period following the train period by:
```{r scorefit}
score_fit(fittmp, !D$trainperiod)$scoreval
score_fit(fittmp, !D$trainperiod)$val
```
In this way it's rather easy to set up different schemes, like optimizing the
......
......@@ -19,7 +19,7 @@ makeit <- function(nam, openit=FALSE, clean=TRUE){
#
unlink(paste0(dirnam,"tmp-setup-data/"), recursive=TRUE)
makeit("setup-data", openit=FALSE)
makeit("setup-data", openit=TRUE)
#
unlink(paste0(dirnam,"tmp-setup-and-use-model/"), recursive=TRUE)
......
......@@ -441,7 +441,7 @@ See the help `?make_tday` for more details.
If we want to use observations in inputs to a model, we can use e.g.:
```{r}
D$Tao <- make_input(D$Ta.obs, kseq=1:36)
D$Tao <- make_input(D$Taobs, kseq=1:36)
model$add_inputs(Tao = "lp(Tao, a1=0.99)")
```
......
......@@ -284,7 +284,7 @@ observations) can be generated by:
i <- 200:296
plot(D$t[i], D$I$k8[i], type="l", col=2, xlab="Time", ylab="Global radiation (W/m²)")
## Add the observations
lines(D$t[i], D$I.obs[i])
lines(D$t[i], D$Iobs[i])
legend("topright", c("8-step forecasts","Observations"), bg="white", lty=1, col=2:1)
```
......@@ -292,7 +292,7 @@ Notice how the are not aligned, since the forecasts are 8 hours ahead. To align
them the forecasts must be lagged 8 steps by:
```{r}
plot(D$t[i], lagvec(D$I$k8[i], 8), type="l", col=2, xlab="Time", ylab="Global radiation (W/m²)")
lines(D$t[i], D$I.obs[i])
lines(D$t[i], D$Iobs[i])
legend("topright", c("8-step forecasts lagged","Observations"), bg="white", lty=1, col=2:1)
```
......@@ -341,7 +341,7 @@ So lagging (thus aligning in time) makes less slightly less scatter.
A wrapper for the `pairs` function is provided for a `data.list`, which can
generate very useful explorative plots:
```{r, fig.height=figwidth}
pairs(D, nms=c("heatload","Ta.obs","Ta","t"), kseq=c(1,8,24))
pairs(D, nms=c("heatload","Taobs","Ta","t"), kseq=c(1,8,24))
```
Note how the sequence of included horizons are specified in the `kseq` argument,
and note that the forecasts are lagged to be aligned in time. See `?pairs.data.list` for more details.
......@@ -354,7 +354,7 @@ x <- lagvec(D$I$k8, 8)
## Take a smaller range
x <- x[i]
## Take the observations
y <- D$I.obs[i]
y <- D$Iobs[i]
## Fit a linear regression model
fit <- lm(y ~ x)
## Plot the result
......@@ -429,6 +429,3 @@ Dsub2 <- as.data.list(Df)
summary(Dsub2)
summary(Dsub)
```
## Literature
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment