Skip to content
Snippets Groups Projects
Commit 228d2856 authored by hgb's avatar hgb
Browse files

initial commit of quantileReg - lot of fixes needed to be done

parent a5df283b
Branches
No related tags found
No related merge requests found
...@@ -15,7 +15,8 @@ Imports: ...@@ -15,7 +15,8 @@ Imports:
R6 (>= 2.2.2), R6 (>= 2.2.2),
splines (>= 3.1.1), splines (>= 3.1.1),
pbs, pbs,
digest digest,
quantreg
LinkingTo: Rcpp, RcppArmadillo LinkingTo: Rcpp, RcppArmadillo
Suggests: Suggests:
knitr, knitr,
......
...@@ -36,3 +36,15 @@ rls_update_cpp <- function(y, X, theta, P, lambda, k, n, np, istart, kmax) { ...@@ -36,3 +36,15 @@ rls_update_cpp <- function(y, X, theta, P, lambda, k, n, np, istart, kmax) {
.Call('_onlineforecast_rls_update_cpp', PACKAGE = 'onlineforecast', y, X, theta, P, lambda, k, n, np, istart, kmax) .Call('_onlineforecast_rls_update_cpp', PACKAGE = 'onlineforecast', y, X, theta, P, lambda, k, n, np, istart, kmax)
} }
rq_init_cpp <- function(X, R, beta, n, weights) {
.Call('_onlineforecast_rq_init_cpp', PACKAGE = 'onlineforecast', X, R, beta, n, weights)
}
rq_simplex_cpp <- function(X, Ih, Ihc, IH, K, n, xB, P, tau) {
.Call('_onlineforecast_rq_simplex_cpp', PACKAGE = 'onlineforecast', X, Ih, Ihc, IH, K, n, xB, P, tau)
}
rq_update_cpp <- function(newX, X, IX, Iy, Ih, Ihc, beta, Rny, K, n, xB, P, tau, k, kk, i, n_in_bin, W = 1) {
.Call('_onlineforecast_rq_update_cpp', PACKAGE = 'onlineforecast', newX, X, IX, Iy, Ih, Ihc, beta, Rny, K, n, xB, P, tau, k, kk, i, n_in_bin, W)
}
...@@ -18,3 +18,26 @@ ...@@ -18,3 +18,26 @@
#' } #' }
#' @source See \url{https://onlineforecasting.org/examples/datasets.html}. #' @source See \url{https://onlineforecasting.org/examples/datasets.html}.
"Dbuilding" "Dbuilding"
#' Observations and weather forecasts from a PV plant, weather station and Danish Meteorological Institute (DMI)
#'
#' Data of the period from 2008-12-29 to 2011-01-01.
#'
#' Hourly average values. The time point is set in the end of the hour.
#'
#' Set in the format of a data.list used as input to forecast models in the onlineforecast package.
#'
#' @format A data list with 1854 rows and 7 variables:
#' \describe{
#' \item{t}{Time in GMT as POSIXct}
#' \item{Ps}{The generated power of a PV plant in MW}
#' \item{Iobs}{Observed global radiation at the weather station in W/m^2}
#' \item{I}{Weather forecasts of global radiation up to 36 hours ahead from DMI in W/m^2}
#' \item{Ta}{Weather forecasts of ambient temperature up to 36 hours ahead from DMI in Celcius}
#' \item{Ws}{Weather forecasts of wind speed up to 36 hours ahead from DMI in m/s}
#' \item{Wd}{Weather forecasts of wind direction up to 36 hours ahead from DMI in degrees}
#' \item{tday}{The time of the day from the given prediction time at the certain prediction horizon in hours}
#' }
#' @source See \url{https://onlineforecasting.org/examples/datasets.html}.
"Dsolar"
\ No newline at end of file
pinball_loss <- function(r, tau){
return(sum(ifelse(r < 0, (tau - 1)*r, tau*r)))
}
qmodel <- R6::R6Class("qmodel", public = list(
####### Store Data and information for algorithm #####
info = list(),
####### Store Coefficients Data #######
beta = list(),
####### Fixed Parameters ######
## Weighted
W = NA,
n_in_bin = NA,
## Number of predictiors
K = NA,
## Quantiles
tau = NA,
## Cold start size
N1 = NA,
debug = NA,
####### Quantile Information ########
## Index of the columns of the X matrix of Xny
IX = NA,
## Index of the column of y
Iy = NA,
## The design matrix
X = list(),
####### debug Information ########
listIH = list(),
####### Forecast Information ########
inputs = list(),
kseq = NA,
output = NA,
prm = NA,
regprmexpr = NA,
regprm = NA,
## Offline parameters
prmbounds = as.matrix(data.frame(lower=NA, init =NA, upper=NA)),
datatr = NA,
qrFIT = NA,
maxlagAR = NA,
yAR = NA,
Ypred = NA,
initialize = function(N1 = NULL, debug = FALSE){
if(is.null(N1)) stop("The number of data points for the cold start needs to be defined")
self$N1 <- N1
self$debug <- debug
},
#----------------------------------------------------------------
# Get the transformation parameters
get_prmbounds = function(nm){
if(nm == "init"){
if(is.null(dim(self$prmbounds))){
val <- self$prmbounds[nm]
}else{
val <- self$prmbounds[ ,nm]
if(is.null(nams(val))){
nams(val) <- rownames(self$prmbounds)
}
}
}
if(nm == "lower"){
if("lower" %in% nams(self$prmbounds)){
val <- self$prmbounds[,"lower"]
if(is.null(nams(val))){
nams(val) <- rownames(self$prmbounds)
}
}else{
val <- -Inf
}
}
if(nm == "upper"){
if("upper" %in% nams(self$prmbounds)){
val <- self$prmbounds[,"upper"]
if(is.null(nams(val))){
nams(val) <- rownames(self$prmbounds)
}
}else{
val <- Inf
}
}
names(val) <- row.names(self$prmbounds)
return(val)
},
#----------------------------------------------------------------
#----------------------------------------------------------------
# Add the transformation parameters and bounds for optimization
add_prmbounds = function(...) {
dots <- list(...)
for (i in 1:length(dots)) {
nm <- names(dots)[i]
if (nm %in% rownames(self$prmbounds)) {
self$prmbounds[nm, ] <- dots[[i]]
} else {
if(nrow(self$prmbounds) == 1 & is.na(self$prmbounds[1,2])){
self$prmbounds[1, ] <- dots[[i]]
}else{
self$prmbounds <- rbind(self$prmbounds, dots[[i]])
}
rownames(self$prmbounds)[nrow(self$prmbounds)] <- nm
}
}
},
#----------------------------------------------------------------
#----------------------------------------------------------------
# Resets the input states
reset_state = function(){
# Reset the inputs state
lapply(self$inputs, function(input){
input$state_reset()
})
# Reset stored data
self$datatr <- NA
#self$yAR <- NA
},
#----------------------------------------------------------------
#----------------------------------------------------------------
# Insert the transformation parameters prm in the input expressions and regression expressions, and keep them (simply string manipulation)
insert_prm = function(prm){
# If just NA or NULL given, then don't do anything
if(is.null(prm) | (is.na(prm)[1] & length(prm) == 1)){
return(NULL)
}
# MUST INCLUDE SOME checks here and print useful messages if something is not right
if(any(is.na(prm))){ stop(pst("None of the parameters (in prm) must be NA: prm=",prm)) }
# Keep the prm
self$prm <- prm
# Find if any opt parameters, first the one with "__" hence for the inputs
pinputs <- prm[grep("__",nams(prm))]
# If none found for inputs, then the rest must be for regression
if (length(pinputs) == 0 & length(prm) > 0) {
preg <- prm
} else {
preg <- prm[-grep("__",nams(prm))]
}
# ################
# For the inputs, insert from prm if any found
if (length(pinputs)) {
pnms <- unlist(getse(strsplit(nams(pinputs),"__"), 1))
pprm <- unlist(getse(strsplit(nams(pinputs),"__"), 2))
#
for(i in 1:length(self$inputs)){
for(ii in 1:length(pnms)){
# Find if the input i have prefix match with the opt. parameter ii
if(pnms[ii]==nams(self$inputs)[i]){
# if the opt. parameter is in the expr, then replace
self$inputs[[i]]$expr <- private$replace_value(name = pprm[ii],
value = pinputs[ii],
expr = self$inputs[[i]]$expr)
}
}
}
}
# ################
# For the fit parameters, insert from prm if any found
if (length(preg) & any(!is.na(self$regprmexpr))) {
nams(preg)
for(i in 1:length(preg)){
# if the opt. parameter is in the expr, then replace
self$regprmexpr <- private$replace_value(name = nams(preg)[i],
value = preg[i],
expr = self$regprmexpr)
}
}
self$regprm <- eval(parse(text = self$regprmexpr))
},
#----------------------------------------------------------------
#----------------------------------------------------------------
# Add inputs to the model
add_inputs = function(...){
dots <- list(...)
for (i in 1:length(dots)){
self$inputs[[ nams(dots)[i] ]] <- onlineforecast:::input_class$new(dots[[i]], model=self)
}
},
#----------------------------------------------------------------
# Add the expression (as character) which generates the regression parameters
add_regprm = function(regprmexpr){
self$regprmexpr <- regprmexpr
self$regprm <- eval(parse(text = self$regprmexpr))
},
#----------------------------------------------------------------
#----------------------------------------------------------------
# Function for transforming the input data to the regression data
transform_data = function(data){
# Evaluate for each input the expresssion to generate the model input data
L <- lapply(self$inputs, function(input){
# Evaluate the expression (input$expr)
L <- input$evaluate(data)
# Must return a list
if(class(L)[1]=="matrix"){ return(list(as.data.frame(L))) }
if(class(L)[1]=="data.frame"){ return(list(L)) }
if(class(L)[1]!="list"){ stop(pst("The value returned from evaluating: ",input$expr,", was not a matrix, data.frame or a list of them."))}
if(class(L[[1]])[1]=="matrix"){ return(lapply(L, function(mat){ return(as.data.frame(mat)) })) }
return(L)
})
# Put together in one data.list
L <- structure(do.call(c, L), class="data.list")
#
return(L)
}
),
# Private functions
private = list(
#----------------------------------------------------------------
#----------------------------------------------------------------
# Replace the value in "name=value" in expr
replace_value = function(name, value, expr){
# First make regex
pattern <- gsub("\\.", ".*", name)
# Try to find it in the input
pos <- regexpr(pattern, expr)
# Only replace if prm was found
if(pos>0){
pos <- c(pos+attr(pos,"match.length"))
# Find the substr to replace with the prm value
(tmp <- substr(expr, pos, nchar(expr)))
pos2 <- regexpr(",|)", tmp)
# Insert the prm value and return
expr <- pst(substr(expr,1,pos-1), "=", value, substr(expr,pos+pos2-1,nchar(expr)))
# Print? Not used now
#if(printout){ message(names(value),"=",value,", ",sep="")}
}
return(expr)
}
#----------------------------------------------------------------
))
#' @title quantile_fit
#'
#' @description This function plots the quantiles
#' @param Lfor List of the data needed to plot
#' @export
#' @examples
#' plotprob()
quantile_fit <- function(prm, model, q, data, printout = TRUE, returnanalysis = TRUE, multistep = NA){
if(printout){
message("----------------")
if(is.na(prm[1])){
message("prm=NA, so current parameters are used.")
}else{
print_to_message(prm)
}
}
# First insert the prm into the model input expressions
model$insert_prm(prm)
model$n_in_bin <- round(2*(1/(1-model$regprm$lambda)))
## Annoying!
if(model$n_in_bin < model$N1){
model$n_in_bin = model$N1 + 1
}
## Transform the data
# Since rls_fit is run from scratch, the init the stored inputs data (only needed when running iteratively)
#model$datatr <- NA
#model$yAR <- NA
# Reset the model state (e.g. inputs state, stored iterative data, ...)
model$reset_state()
if(q %in% as.numeric(sub(".*q", "", names(model$info)))){
model$info <- list()
model$beta <- list()
}
# Generate the 2nd stage inputs (i.e. the transformed data)
model$datatr <- model$transform_data(data)
## Fit the cold start
#X <- data.frame(matrix(unlist(model$datatr), ncol = length(model$datatr), byrow = F))
#colnames(X) <- names(model$datatr)
Ypred <- NULL
for(k in model$kseq){
cat("horizon :", k, "\n")
cat("n_in_bin :", model$n_in_bin, "\n")
## Setup the design matrix for k step
X <- as.data.frame(subset(model$datatr, kseq=k))
## Lag them to match to the
X <- onlineforecast:::lagdf.matrix(X, k)
X$y <- data[[model$output]]
X <- X[-1:-k,]
fit <- rq(paste("y ~ -1 + ", paste(names(X[, -which(names(X) %in% c("y"))]), collapse = " + "), sep = ""), tau = q, data = X[1:(model$N1-k),])
beta <- as.matrix(fit$coefficients)
r <- as.matrix(fit$residuals, ncol = length(q))
for(i in 1:dim(r)[2]){
## So these guys are the zeros, or lie on the constraints?
IdxNull <- which(abs(r[,i]) < 10^-2)
r[IdxNull, i] <- 0
}
### Initialization of algorithm
model$info[[paste0("q", q)]][[paste0("k", k)]] <- rq_init_cpp(X = as.matrix(X[1:(model$N1-k),]),
R = r,
beta = beta,
n = model$N1-k,
weights = model$regprm$lambda^rev(seq(1,length(r),1)))
## The quantile estimation
Ypredres <- update_rq(as.matrix(X[(model$N1+1 - k):nrow(X),]), q, k, model = model)
### Get the results, Ypred, loss
Ypred <- cbind(Ypred, Ypredres)
}
## Compute the loss function, maybe the residuals are compute in different place
loss <- sapply(model$kseq, function(i) pinball_loss(r = model$info[[paste0("q", q)]][[paste0("k", i)]]$R[data$scoreperiod[(i+1):length(data$scoreperiod)]], q))
model$Ypred <- Ypred
if(!returnanalysis){
#Yhat <- quantile_prediction(model = model, multistep = model$kseq)
cat("Quantile, tau:", q, "\n")
cat("the loss: ", loss, "\n")
val <- sum(loss)
if(printout){
print_to_message(c(loss,sum=val))
}
return(val)
} else{
val <- sum(loss)
if(printout){
print_to_message(c(loss,sum=val))
}
## Get predictions
#Yhat <- model$getPred(multistep = multistep)
#Yhat <- quantile_prediction(model = model, multistep = multistep)
#invisible(structure(list(Yhat = Yhat, model = model, data = data, scoreval = loss), class = c("quantile_fit")))
}
}
print_to_message <- function(...) {
message(paste(utils::capture.output(print(...)), collapse="\n"))
}
#' @title quantile_optim
#'
#' @description This function plots the quantiles
#' @param Lfor List of the data needed to plot
#' @export
#' @examples
#' plotprob()
quantile_optim <- function(model, data){
# Take the parameters bounds from the parameter bounds set in the model
init <- model$get_prmbounds("init")
lower <- model$get_prmbounds("lower")
upper <- model$get_prmbounds("upper")
# If bounds are NA, then set
if(any(is.na(lower))){ lower[is.na(lower)] <- -Inf}
if(any(is.na(upper))){ lower[is.na(upper)] <- Inf}
resList <- list()
for(i in (1:length(model$tau))){
cat("Quantile, tau: ", model$tau[i], "\n")
resList[[i]] <- optim(par = init,
fn = quantile_fit,
# Parameters to pass to quantile_fit
model = model,
q = model$tau[i],
data = Dnew,
returnanalysis = FALSE,
# Parameters to pass to optim
lower = lower,
upper = upper,
method = "L-BFGS-B")
}
#model$optPar <- sapply(1:length(model$tau), function(q) return(resList[[i]]$par))
return(resList)
}
#' @title quantile_prediction
#'
#' @description This function plots the quantiles
#' @param Lfor List of the data needed to plot
#' @export
#' @examples
#' plotprob()
quantile_predict <- function(model, datatr){
Yhat <- lapply(model$tau, function(q){
yHat <- sapply(model$kseq, function(k){
## Setup the design matrix for k step
X <- as.data.frame(subset(datatr, kseq=k))
## Lag them to match to the
#X <- onlineforecast:::lagdf.matrix(X, k)
yhat <- as.numeric(rep(NA, nrow(X)))
for(i in ((k):nrow(X))){
if(i <= (model$N1)) {
j <- 1
} else{
j <- i - model$N1
}
#browser()
yhat[i] <- t(as.numeric(X[i,])) %*% model$beta[[paste0("q",q)]][[paste0("k",k)]][j,]
}
return(yhat)
})
nams(yHat) <- pst("k", model$kseq)
return(yHat)
})
names(Yhat) <- pst("q", model$tau)
return(Yhat)
}
quantile_update <- function(model, datatr = NA, y = NA){
# Find the number of parameters for the regression
np <- length(datatr)
# Check if data was kept
kept_input_data <- !is.na(model$datatr[1])
# Check if data was kept
kept_input_data <- !is.na(model$datatr[1])
# Parameters for quantile
lambda <- model$regprm$lambda
## n_in_bin could be change!!
n_in_bin <- round(2*(1/(1-lambda)))
model$run(X = X[, -which(names(X) %in% c("y"))],
Y = X[, "y"],
tau = tau,
residuals = as.matrix(r),
beta = matrix(beta, nrow = length(datatr)),
Time = data[["t"]],
n_in_bin = round(2*(1/(1-model$regprm$lambda))),
Weights = model$regprm$lambda)
invisible()
}
#' @title quantile_fit
#'
#' @description This function plots the quantiles
#' @param Lfor List of the data needed to plot
#' @export
#' @examples
#' plotprob()
sortIdx = function(val){
newVal <- sort(val)
idxVal <- order(val)
return(list(newVal = newVal, idxVal = idxVal))
}
#' @title update_rq
#'
#' @description This function plots the quantiles
#' @param Lfor List of the data needed to plot
#' @export
#' @examples
#' plotprob()
update_rq <- function(Xny, tau, k, model = NULL, debug = TRUE){
if(is.null(model)) stop("Model needs to be initialized before update")
N <- nrow(Xny)
res <- model$info[[paste0("q", tau)]][[paste0("k", k)]]
## Maybe this is no correct - we should put weights on this!
Xold <- res$X
# Prepare for keeping for the parameter estimates
BETA <- matrix(as.numeric(NA), nrow = N, ncol = model$K)
i <- 1
counter <- 1
Ypred <- matrix(as.numeric(NA), nrow = N, ncol = 1)
if(debug) {
print("yeah")
list_idx <- matrix(as.numeric(NA), nrow = 1, ncol = length(res$Ih)+1)
list_idx[1,] <- c(0,res$Ih)
}
while(counter <= N){
#cat("counter", counter, "\n")
beta <- res$xB[1:model$K]
BETA[counter,] <- t(beta)
j <- 0
res <- rq_update_cpp(newX = matrix(Xny[counter,], nrow = 1),
X = Xold,
IX = model$IX,
Iy = model$Iy,
Ih = res$Ih,
Ihc = res$Ihc,
beta = beta,
Rny = res$R,
K = model$K,
n = nrow(Xold),
xB = res$xB,
P = res$P,
tau = tau,
k = 0,
kk = 0,
i = i,
n_in_bin = model$n_in_bin,
W = model$regprm$lambda)
Ypred[counter] <- res$Ypred
i = res$i
Xold <- res$X
res$IH <- res$Ih
#cat("Here", "\n")
#if(counter == 119) browser()
if(debug) {
list_idx <- rbind(list_idx, c(counter, res$Ih))
}
resAlg <- rq_simplex_cpp(X = Xold[,1:model$K],
Ih = res$Ih,
Ihc = res$Ihc,
IH = as.matrix(res$Ih),
K = model$K,
n = nrow(Xold),
xB = res$xB,
P = res$P,
tau = tau)
res$CON <- resAlg$CON
res$s <- resAlg$s
res$g <- resAlg$g
res$gain <- resAlg$gain
res$md <- resAlg$md
res$alpha <- resAlg$alpha
res$h <- resAlg$h
res$IH <- resAlg$IH
res$cq <- resAlg$cq
res$q <- resAlg$q
if(length(res$md) == 0) res$md <- 1
while(res$gain <= 0 & res$md < 0 & j < 24 & res$CON < 10^6){
j <- j + 1
z = res$xB - as.numeric(res$alpha) * res$h
IhM <- res$Ih[res$s + 1]
IhcM <- res$Ihc[res$q + 1]
res$Ih[res$s + 1] <- IhcM
res$Ihc[res$q+1] <- IhM
res$P[res$q+1] <- res$cq
res$xB <- z
res$xB[res$q + 1 + model$K] <- res$alpha
dummy <- sortIdx(res$Ih)
res$Ih <- dummy$newVal; IndexIh <- dummy$idxVal
dummy <- sortIdx(res$Ihc)
res$Ihc <- dummy$newVal; IndexIhc <- dummy$idxVal
if(debug) {
list_idx <- rbind(list_idx, c(counter, res$Ih))
}
res$P <- res$P[IndexIhc]
xBm <- res$xB[((model$K+1):length(res$xB))]
xBm <- xBm[IndexIhc]
res$xB[((model$K+1):length(res$xB))] <- xBm
#cat("Here2", "\n")
resAlg <- rq_simplex_cpp(X = Xold[,1:model$K],
Ih = res$Ih,
Ihc = res$Ihc,
IH = as.matrix(res$Ih),
K = model$K,
n = nrow(Xold),
xB = res$xB,
P = res$P,
tau = tau)
res$CON <- resAlg$CON
res$s <- resAlg$s
res$g <- resAlg$g
res$gain <- resAlg$gain
res$md <- resAlg$md
res$alpha <- resAlg$alpha
res$h <- resAlg$h
res$IH <- resAlg$IH
res$cq <- resAlg$cq
res$q <- resAlg$q
if(length(res$md) == 0) res$md <- 1
}
counter = counter + 1
}
model$info[[paste0("q", tau)]][[paste0("k", k)]] <- res
model$beta[[paste0("q", tau)]][[paste0("k", k)]] <- rbind(model$beta[[paste0("q", tau)]][[paste0("k", k)]], BETA)
if(debug) model$listIH[[paste0("q", tau)]][[paste0("k", k)]] <- list_idx
return(Ypred)
}
File added
### Clean up the workspace
rm(list = ls())
## additional libraries
library(devtools)
library(splines)
library(quantreg)
## Load the onlineforecast package - this is the development version
load_all(".") # need to be under the package directory
# Load the data
D <- Dsolar
names(D)
# Short Data Anlysis
tstart <- as.POSIXct("2009-01-01", tz = "GMT")
tend <- as.POSIXct("2011-01-01", tz = "GMT")
D <- subset(D, in_range(tstart,D$t,tend))
plotmat <- matrix(1:3, nrow = 3, byrow=TRUE)
layout(plotmat, widths=c(1,1,1), heights=c(1,1,1))
par(mar=c(0, 3.5, 0.5, 1), lwd = 2, cex.lab = 1.5, cex.axis = 1.5)
plot(D$t, D$Ps, type = "l", bty = "l", xlab = "", ylab = "", axes = F)
axis(2, seq(0,ceiling(max(D$Ps)+0.2), by = 0.5), las = 1)
plot(D$t, lagvec(D$I$k1,1), type = "l", col = "gold", bty = "l", xlab = "", ylab = "", axes = F)
lines(D$t, D$Iobs, col = rgb(0, 0, 255, max = 255, alpha = 125))
axis(2, seq(0,ceiling(max(D$I$k1)+0.2), by = 100), las = 1)
plot(D$t, lagvec(D$Ta$k1,1), type = "l", col = rgb(240, 37, 37, max = 255, alpha = 150), bty = "l", xlab = "", ylab = "", axes = F, ylim = c(min(D$Ta$k1) - 4,max(D$Ta$k1)))
axis(2, seq(floor(min(D$Ta$k1)),ceiling(max(D$Ta$k1)+0.2), by = 10), las = 1)
axis.POSIXct(side = 1, x = D$t, xaxt = "s",
at = seq(D$t[1], D$t[length(D$t)], by = "1 month"),
format = "%Y-%m-%d", mgp = c(4,-1,-2))
# Make a training set of 3 months, and a then test set with 1 month
tstart <- "2009-06-01"
ttest <- "2009-03-31"
tend <- "2009-07-01"
# Cut only the necessary period
D$scoreperiod <- in_range("2009-06-10", D$t)
Dtrain <- subset(D, in_range(tstart,D$t,tend))
model <- qmodel$new(N1 = 50)
model$output = "Ps"
model$add_inputs(I = "bspline(tday, df=5, Boundary.knots=c(5,18)) %**% I")
model$add_regprm("rls_prm(lambda=0.9)")
model$add_prmbounds(lambda = c(0.99, 0.999, 0.9999))
## First 5 columns are the input Matrix
model$IX <- 0:4
## The output column is the 6th one
model$Iy <- 5
model$K <- 5
model$datatr <- model$transform_data(Dtrain)
model$tau <- c(0.05,0.2, 0.5, 0.8,0.95)
model$kseq <- 1:24
#opt_model <- quantile_optim(model = model, data = Dnew)
plot(Dtrain$I$k1[1:100], type = "l")
lines(model$datatr$I.bs1$k1, col = "red")
lines(model$datatr$I.bs2$k1, col = "red")
lines(model$datatr$I.bs3$k1, col = "red")
lines(model$datatr$I.bs4$k1, col = "red")
lines(model$datatr$I.bs5$k1, col = "red")
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[1])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[2])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[3])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[4])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[5])
Pred_model <- quantile_predict(model = model, datatr = model$datatr)
str(Pred_model)
k <- 24
colSeq <- grey(seq(0.9,0.1,len=3))
#colSeq <- colorJet(nPolygons)
require( RColorBrewer )
colSeq <- brewer.pal(11 , "Spectral" )
for(i in k:(length(Dtrain$t)-1)){
par(mfrow = c(2,1))
par(mar = c(2,4,2,2)) # bottom left top right
with(Pred_model, {
plot(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], bty = "l", lwd = 2, col = "black", pch=19, cex=0.5, axes = FALSE, xaxt = "n",
type="n", ylim = range(Dtrain$Ps[(i-10):(i+k)], q0.5[i,], q0.95[i,], q0.05[i,], na.rm = T), main = "TimeAdaptive with Weights", xlab = "Time", ylab = "Ps")
axis(2)
polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
c(Pred_model$q0.05[i,],rev(Pred_model$q0.2[i,])), col= colSeq[3], border = NA)
polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
c(Pred_model$q0.8[i,],rev(Pred_model$q0.95[i,])), col= colSeq[3], border = NA)
polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
c(Pred_model$q0.2[i,],rev(Pred_model$q0.5[i,])), col= colSeq[4], border = NA)
polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
c(Pred_model$q0.5[i,],rev(Pred_model$q0.8[i,])), col= colSeq[4], border = NA)
lines(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], col = "black", type = "b", lwd = 2)
lines(Dtrain$t[(i+1):(i+k)], q0.5[i,], type = "b", col = "grey", lwd = 2)
axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)],
at = seq(from = Dtrain$t[(i-10):(i+k)][1],
to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
las = 1, cex.axis = 1, srt = 45)
#lines(Pred_model$q0.05[i,], type = "b", col = "blue")
#lines(Pred_model$q0.95[i,], type = "b", col = "blue")
plot(Dtrain$t[(i-10):(i+k)], c(rep(NA,11), as.numeric(Dtrain$I[i, 1:k])), type = "b", col = "steelblue", axes = FALSE, xlab = "Time", ylab = "Temp", lwd = 2, ylim = range(Dtrain$I[i, 1:k], Dtrain$Iobs[(i+1):(i+k)]))
#lines(Dtrain$t[i], Dtrain$Iobs[i+1,1], type = "b", col = "red", axes = FALSE)
lines(Dtrain$t[(i+1):(i+k)], Dtrain$Iobs[(i+1):(i+k)], type = "b", col = "red", lwd = 2)
axis(2)
abline(h=0,v=Dtrain$t[i],lty=2, col = "lightgrey", lwd = 2)
legend("topleft", legend=c("I Obs", "I Pred"),
col=c("red", "steelblue"), lty=1, cex=1, lwd = 2)
i
axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)],
at = seq(from = Dtrain$t[(i-10):(i+k)][1],
to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
las = 1, cex.axis = 1, srt = 45)
})
}
...@@ -43,10 +43,75 @@ BEGIN_RCPP ...@@ -43,10 +43,75 @@ BEGIN_RCPP
return rcpp_result_gen; return rcpp_result_gen;
END_RCPP END_RCPP
} }
// rq_init_cpp
Rcpp::List rq_init_cpp(arma::mat X, arma::vec R, arma::vec beta, unsigned int n, arma::vec weights);
RcppExport SEXP _onlineforecast_rq_init_cpp(SEXP XSEXP, SEXP RSEXP, SEXP betaSEXP, SEXP nSEXP, SEXP weightsSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP);
Rcpp::traits::input_parameter< arma::vec >::type R(RSEXP);
Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP);
Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
Rcpp::traits::input_parameter< arma::vec >::type weights(weightsSEXP);
rcpp_result_gen = Rcpp::wrap(rq_init_cpp(X, R, beta, n, weights));
return rcpp_result_gen;
END_RCPP
}
// rq_simplex_cpp
Rcpp::List rq_simplex_cpp(arma::mat X, arma::uvec Ih, arma::uvec Ihc, arma::mat IH, unsigned int K, unsigned int n, arma::vec xB, arma::vec P, double tau);
RcppExport SEXP _onlineforecast_rq_simplex_cpp(SEXP XSEXP, SEXP IhSEXP, SEXP IhcSEXP, SEXP IHSEXP, SEXP KSEXP, SEXP nSEXP, SEXP xBSEXP, SEXP PSEXP, SEXP tauSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP);
Rcpp::traits::input_parameter< arma::uvec >::type Ih(IhSEXP);
Rcpp::traits::input_parameter< arma::uvec >::type Ihc(IhcSEXP);
Rcpp::traits::input_parameter< arma::mat >::type IH(IHSEXP);
Rcpp::traits::input_parameter< unsigned int >::type K(KSEXP);
Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
Rcpp::traits::input_parameter< arma::vec >::type xB(xBSEXP);
Rcpp::traits::input_parameter< arma::vec >::type P(PSEXP);
Rcpp::traits::input_parameter< double >::type tau(tauSEXP);
rcpp_result_gen = Rcpp::wrap(rq_simplex_cpp(X, Ih, Ihc, IH, K, n, xB, P, tau));
return rcpp_result_gen;
END_RCPP
}
// rq_update_cpp
Rcpp::List rq_update_cpp(arma::mat newX, arma::mat X, arma::uvec IX, arma::uvec Iy, arma::uvec Ih, arma::uvec Ihc, arma::vec beta, arma::vec Rny, unsigned int K, unsigned int n, arma::vec xB, arma::vec P, double tau, unsigned int k, arma::uvec kk, unsigned int i, unsigned int n_in_bin, double W);
RcppExport SEXP _onlineforecast_rq_update_cpp(SEXP newXSEXP, SEXP XSEXP, SEXP IXSEXP, SEXP IySEXP, SEXP IhSEXP, SEXP IhcSEXP, SEXP betaSEXP, SEXP RnySEXP, SEXP KSEXP, SEXP nSEXP, SEXP xBSEXP, SEXP PSEXP, SEXP tauSEXP, SEXP kSEXP, SEXP kkSEXP, SEXP iSEXP, SEXP n_in_binSEXP, SEXP WSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::mat >::type newX(newXSEXP);
Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP);
Rcpp::traits::input_parameter< arma::uvec >::type IX(IXSEXP);
Rcpp::traits::input_parameter< arma::uvec >::type Iy(IySEXP);
Rcpp::traits::input_parameter< arma::uvec >::type Ih(IhSEXP);
Rcpp::traits::input_parameter< arma::uvec >::type Ihc(IhcSEXP);
Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP);
Rcpp::traits::input_parameter< arma::vec >::type Rny(RnySEXP);
Rcpp::traits::input_parameter< unsigned int >::type K(KSEXP);
Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
Rcpp::traits::input_parameter< arma::vec >::type xB(xBSEXP);
Rcpp::traits::input_parameter< arma::vec >::type P(PSEXP);
Rcpp::traits::input_parameter< double >::type tau(tauSEXP);
Rcpp::traits::input_parameter< unsigned int >::type k(kSEXP);
Rcpp::traits::input_parameter< arma::uvec >::type kk(kkSEXP);
Rcpp::traits::input_parameter< unsigned int >::type i(iSEXP);
Rcpp::traits::input_parameter< unsigned int >::type n_in_bin(n_in_binSEXP);
Rcpp::traits::input_parameter< double >::type W(WSEXP);
rcpp_result_gen = Rcpp::wrap(rq_update_cpp(newX, X, IX, Iy, Ih, Ihc, beta, Rny, K, n, xB, P, tau, k, kk, i, n_in_bin, W));
return rcpp_result_gen;
END_RCPP
}
static const R_CallMethodDef CallEntries[] = { static const R_CallMethodDef CallEntries[] = {
{"_onlineforecast_lp_vector_cpp", (DL_FUNC) &_onlineforecast_lp_vector_cpp, 2}, {"_onlineforecast_lp_vector_cpp", (DL_FUNC) &_onlineforecast_lp_vector_cpp, 2},
{"_onlineforecast_rls_update_cpp", (DL_FUNC) &_onlineforecast_rls_update_cpp, 10}, {"_onlineforecast_rls_update_cpp", (DL_FUNC) &_onlineforecast_rls_update_cpp, 10},
{"_onlineforecast_rq_init_cpp", (DL_FUNC) &_onlineforecast_rq_init_cpp, 5},
{"_onlineforecast_rq_simplex_cpp", (DL_FUNC) &_onlineforecast_rq_simplex_cpp, 9},
{"_onlineforecast_rq_update_cpp", (DL_FUNC) &_onlineforecast_rq_update_cpp, 18},
{NULL, NULL, 0} {NULL, NULL, 0}
}; };
......
#include "RcppArmadillo.h"
#include <Rcpp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
Rcpp::List rq_init_cpp(arma::mat X,
arma::vec R,
arma::vec beta,
unsigned int n,
arma::vec weights){
//arma::vec idx(n);
arma::mat L, U, Pi;
unsigned int Lr0 = count(R.begin(), R.end(), 0);
//unsigned int idx;
if(Lr0 > beta.n_elem){
Rcout << "Worked" << "\n";
arma::uvec Irs = sort_index(abs(R));
Irs = Irs.subvec(0,Lr0-1);
arma::mat Xh = X.rows(Irs);
lu(L, U, Pi, Xh);
arma::vec In = arma::linspace(0,Lr0,Lr0+1);
arma::vec rI = arma::linspace(0, (Lr0)-beta.n_elem-1, (Lr0)-beta.n_elem);
for(unsigned int i = (beta.n_elem); i < Lr0; i++){
//idx = In(find(Pi.row(i) == 1));
R(Irs(find(Pi.row(i) == 1))).fill(1.e-15);
}
//R(Irs(rI)) = 1.e-15;
}
// Find all index where r == 0 - Vertices
arma::uvec Ih = find(R == 0);
arma::uvec Ihc = find(abs(R) > 0);
arma::vec P = sign(R.elem(Ihc));
R.elem(find(abs(R) < 1.e-15)).fill(0);
arma::vec W = weights % R;
arma::vec xB = join_cols(beta, abs(W(Ihc)));
return List::create(Named("xB") = xB,
Named("Ih") = Ih,
Named("Ihc") = Ihc,
Named("P") = P,
Named("R") = R,
Named("X") = X);
}
\ No newline at end of file
#include "RcppArmadillo.h"
#include <Rcpp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
Rcpp::List rq_simplex_cpp(arma::mat X,
arma::uvec Ih,
arma::uvec Ihc,
arma::mat IH,
unsigned int K,
unsigned int n,
arma::vec xB,
arma::vec P,
double tau){
arma::mat invXh = inv(X.rows(Ih));
arma::vec cB = arma::conv_to<arma::vec>::from(P < 0) + P * tau;
arma::vec cC = arma::join_cols(arma::ones<arma::vec>(K) * tau, arma::ones<arma::vec>(K) * (1-tau));
//Rcout << Xny .rows(Ihc) << "\n";
arma::mat IB2 = -(P * arma::ones<arma::vec>(K).t() % X.rows(Ihc))*invXh;
arma::vec g = IB2.t() * cB;
arma::vec d = cC - arma::join_cols(g,-g);
d.elem(find(abs(d) < 1.e-15)).fill(0);
arma::uvec s = sort_index(d);
arma::vec md = sort(d);
s = s.elem(find(md < 0));
md = md.elem(find(md < 0));
arma::vec c = arma::ones<arma::vec>(s.size());
for(unsigned int i = 0; i < s.size(); i++){
if(s(i) >= K){
s(i) = s(i) - K;
c(i) = -1;
}
}
//c.elem(s>K).fill(-1);
arma::mat C(c.size(),c.size(),arma::fill::eye);
C.diag() = c;
arma::mat h = arma::join_cols(invXh.cols(s), IB2.cols(s)) * C;
arma::vec xm = xB.elem(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
xm.elem(find(xm < 0)).fill(0);
arma::mat hm = h.rows(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
arma::vec alpha(s.size());
arma::uvec q(s.size());
arma::vec cq(s.size());
arma::uvec idx2 = arma::linspace<arma::uvec>(0,s.size()-1,s.size());
for(unsigned int k = 0; k < s.size(); k++){
arma::vec sigma = xm;
arma::uvec idx = find(hm.col(k) > 1.e-12);
sigma.elem(idx) = xm.elem(idx) / hm(idx, find(idx2 == k));
sigma(find(hm.col(k) <= 1.e-12)).fill(1.e30);
alpha.elem(find(idx2 == k)).fill(sigma.min());
q.elem(find(idx2 == k)).fill(sigma.index_min());
cq.elem(find(idx2 == k)) = c.elem(find(idx2 == k));
}
arma::vec gain = md % alpha;
arma::vec Mgain = sort(gain);
arma::uvec IMgain = sort_index(gain);
arma::uvec IhMid;
double CON = 1.e30;
unsigned int j = 0;
if(gain.size() < 1){
gain.resize(1);
gain(0) = 1;
} else {
while((CON > 1.e6) & (j < s.size())){
IhMid = Ih;
IhMid(s(IMgain(j))) = Ihc(q(IMgain(j)));
IhMid = sort(IhMid);
if(sum(abs(IH.t() - IhMid.t())) == 0){
CON = 1.e30;
} else{
CON = cond(X.rows(IhMid));
}
s = s(IMgain(j));
q = q(IMgain(j));
cq = cq(IMgain(j));
alpha = alpha(IMgain(j));
IH = arma::join_rows(IH, arma::conv_to<arma::vec>::from(IhMid));
h = h.col(IMgain(j));
gain = gain(IMgain(j));
md = md(IMgain(j));
j += 1;
}
}
return List::create(Named("CON") = CON,
Named("s") = s,
Named("g") = g,
Named("q") = q,
Named("gain") = gain,
Named("md") = md,
Named("alpha") = alpha,
Named("h") = h,
Named("IH") = IH,
Named("cq") = cq);
}
\ No newline at end of file
#include "RcppArmadillo.h"
#include <Rcpp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
Rcpp::List rq_update_cpp(arma::mat newX,
arma::mat X,
arma::uvec IX,
arma::uvec Iy,
arma::uvec Ih,
arma::uvec Ihc,
arma::vec beta,
arma::vec Rny,
unsigned int K,
unsigned int n,
arma::vec xB,
arma::vec P,
double tau,
unsigned int k,
arma::uvec kk,
unsigned int i,
unsigned int n_in_bin,
double W = 1){
// Init
arma::vec xBm;
// If weights are given - then weight the past Xny matrix
// Notice, tho the newest observaiton is not weight - maybe need a fix?
if(W != 1){
X = W * X;
}
X = join_cols(X, newX);
// Rolling window of leaving the design matrix
// i.e. the last observation leaves the matrix
// See Møller 2008 if something else is needed
unsigned int Leav = 0;
arma::uvec sortIdx = sort_index(Ih);
//Rcout << "min" << sortIdx << "\n";
//Rcout << "minV" << sortIdx(0) << "\n";
if(Ih(sortIdx(0)) == Leav){
//Rcout << "Here but still hmm" << "\n";
//Rcout << "Ih" << Ih << "\n";
// Compute the inverse of X(h)
arma::mat invXh = inv(X(Ih, IX));
//Rcout << invXh.col(sortIdx(0)) << "\n";
// Computing the c
arma::vec cB = arma::conv_to<arma::vec>::from(P < 0) + P * tau;
// Computing inverse (B^T) * c_B - P * X(bar(h)) * X(b)^-1 * inverse(X(bar(h)))
double g = arma::as_scalar(-(P % (X(Ihc,IX) * invXh.col(sortIdx(0)))).t() * cB);
int signG = (g > 0) - (g < 0);
// h is the direction or the d in the paper
arma::vec h = join_cols(invXh.col(sortIdx(0)), -(P % (X(Ihc,IX) * invXh.col(sortIdx(0))))) * signG;
// Compute the magitude of the step, alpha
arma::vec sigma = arma::zeros(n - K);
arma::vec hm = h.elem(arma::linspace<arma::uvec>(K,h.size() - 1,h.size() - K));
// The absolute of the residuals which are not equal to zero
arma::vec xBm = xB.elem(arma::linspace<arma::uvec>(K,h.size()-1,h.size()-K));
xBm.elem(find(xBm < 0)).fill(0);
//Rcout << "Here but still hmm" << "\n";
//Rcout << "hm" << hm << "\n";
//Rcout << "Sigma" << sigma << "\n";
//Rcout << "hm" << hm.size() << "\n";
//Rcout << "Sigma" << sigma.size() << "\n";
sigma.elem(find(hm > 1.e-10)) = xBm.elem(find(hm>1.e-10)) / hm.elem(find(hm>1.e-10));
sigma.elem(find(hm<= 1.e-10)).fill(1.e30);
//Rcout << "min idx " << index_min(sigma) << "\n";
unsigned int q = index_min(sigma); //sort_index(sigma).elem(0);
double alpha = arma::as_scalar(sigma.row(q));
arma::vec z = xB - alpha * h;
arma::uvec Ihm = Ih.row(sortIdx(0));
Ih.row(sortIdx(0)) = Ihc.row(q);
Ihc(q) = Ihm(0);
xB = z;
xB(q+K) = alpha;
P(q) = signG + (g == 0);
Ih = sort(Ih);
arma::uvec IndexIhc = sort_index(Ihc);
Ihc = sort(Ihc);
P = P(IndexIhc);
xBm = xB(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
xBm = xBm(IndexIhc);
xB(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K)) = xBm;
beta = xB(arma::linspace<arma::uvec>(0,K-1,K));
}
//Rcout << "IX " << IX << "\n";
//Rcout << "X " << X(kk, IX) << "\n";
//Rcout << "Pred " << arma::as_scalar(X(kk, IX) * beta) << "\n";
double Ypred = arma::as_scalar(newX(kk, IX) * beta);
double rny = arma::as_scalar(newX(kk, Iy)) - Ypred;
//Rcout << "rny" << rny << "\n";
int szR = Rny.size();
Rny.resize(szR+1);
Rny(szR) = rny;
int szP = P.size();
P.resize(szP+1);
int szB = xB.size();
xB.resize(szB+1);
if(rny < 0){
P(szP) = -1;
xB(szB) = -rny;
} else {
P(szP) = 1;
xB(szB) = rny;
}
int sz = Ihc.size();
Ihc.resize(sz+1);
Ihc(sz) = n;
//Rcout << "Ihc" << Ihc << "\n";
if(X.n_rows == n_in_bin){
i += 1;
arma::uvec stay = arma::linspace<arma::uvec>(1,P.size()-1,P.size()-1);
P = P(stay);
Ihc = Ihc(stay);
//arma::vec Pnew = P(stay);
//arma::uvec Ihcnew = Ihc(stay);
X = X.rows(sort(join_cols(Ihc,Ih)));
xBm = xB.elem(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
xBm = xBm(stay);
xB = xB.elem(arma::linspace<arma::uvec>(0,xB.size()-2,xB.size()-1));
xB.elem(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K)) = xBm;
//Ihc = Ihc(stay);
Ihc = Ihc - 1;
Ih = Ih - 1;
} else{
n += 1;
}
return List::create(Named("Ih") = Ih,
Named("Ihc") = Ihc,
Named("xB") = xB,
Named("X") = X,
Named("R") = Rny,
Named("P") = P,
Named("n") = n,
Named("i") = i,
Named("Ypred") = Ypred);
}
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment