R/FcExtremeLearningMachines.R

In Mthrun/TSAT: Time Series Analysis Tools (TSAT)

# res = FcExtremeLearningMachines(DataVec, Time)
#
# DESCRIPTION
# Special case for a multilayer perceptrone feed forward network with backpropagation used for forecasting [Huang et al., 2005] 
#
# INPUT
# DataVec             [1:n] numerical vector of time series data
# Time                [1:n] character vector of Time in the length of data
# SplitAt             Index of row where the DataVec is divided into test and train data. If not given n is used
# ForecastHorizon     Scalar defining the timesteps to forecast ahead
# Predictors          [1:n,1:d] data frame or matrix of d predictores with n values each
# OPTIONAL
# No_HiddenLayers     Number of hidden layers, see nnfor::elm
# Scaled              TRUE: automatic scaling
# No_TrainingNetworks Number of Training Networks, see nnfor::elm
# PlotIt              FALSE (default), do nothing. TRUE: plots the forecast versus the validation set.
# ...                 Further arguments passed on to \code{\link[nnfor]{elm}}
#
# OUTPUT
# Model               Model paramters, see example
# Forecast            [1:n-ForecastHorizon] Forecast generated by the model 
# TestData            [(k+1):n] vector, the part of Response not used in the model
# TestTime            [(k+1):n] vector, time of response not used in the model
# TrainData           [1:k] vector, the part of Response used in the model
# TrainTime           [1:k] vector, time of Training data if given
# TrainingForecast    [1:k] vector, forecasted value using TrainData
#
# DETAILS
# The parameters of hidden nodes (not just the weights connecting inputs to hidden nodes) need not be tuned. These hidden nodes 
# can be randomly assigned and never updated (i.e. they are random projection but with nonlinear transforms), or can be inherited 
# from their ancestors without being changed. In most cases, the output weights of hidden nodes are usually learned in a single step. 
# According to their creators, these models are able to produce good generalization performance and learn thousands of times faster 
# than networks trained using backpropagation [Huang et al., 2005].
#
# NOTE: A Forecast Horizon beyond the test data is only possible if no predictors are given.
#
# Author: MCT

FcExtremeLearningMachines=function(DataVec,Time,SplitAt,Predictors,ForecastHorizon,No_HiddenLayers=NULL,Scaled=TRUE,No_TrainingNetworks=10,PlotEvaluation=FALSE,PlotIt=FALSE,tz="UTC",...){
  

  if(missing(Predictors))
    Predictors=NULL
  else {
    if(!is.null(Predictors)) {
      if(!is.data.frame(Predictors)) {
        if(is.matrix(Predictors) && is.numeric(Predictores)) {
          Predictors = as.data.frame(Predictors)
        } else {
          warning("Predictors should be a data frame or numeric matrix, trying to parse predictors")
          Predictors = as.data.frame(Predictors)
        }
      }
    }
  }
  
  if(missing(SplitAt)) {
    warning('Input for SplitAt was not given. Setting SplitAt to length(DataVec)')
    SplitAt = length(DataVec)
  }
  if(missing(ForecastHorizon)) {
    warning('Input for ForecastHorizon was not given. Setting ForecastHorizon to 1')
    ForecastHorizon = 1
  }
  if(missing(Time)) {
    Time = 1:length(DataVec)
    warning('No input for Time was given, will use 1:length(DataVec) for Time')
  }
  
  inputs = checkInputForecasting(DataVec, Time, SplitAt, ForecastHorizon, PlotIt)
  DataVec = inputs$DataVec
  Time = inputs$Time
  SplitAt = inputs$SplitAt
  ForecastHorizon = inputs$ForecastHorizon
  PlotIt = inputs$PlotIt
  
  if(!is.logical(PlotEvaluation)) {
    warning('Input for PlotEvaluation is not logical. Setting PlotEvaluation to FALSE')
    PlotEvaluation = FALSE
  }
  
  if(!is.logical(Scaled)) {
    warning('Input for Scaled is not logical. Setting Scaled to TRUE')
    Scaled = TRUE
  }
  
  if(isTRUE(PlotEvaluation) && isTRUE(PlotIt)) {
    warning('Can either plot model evaluation or forecast evaluation. Will plot PlotIt (forecast evaluation)')
    PlotEvaluation = FALSE
  }
  
  if(!(is.numeric(No_TrainingNetworks) && length(No_TrainingNetworks) == 1 && No_TrainingNetworks >= 0)) {
    warning('No_TrainingNetworks should be a positive scalar. Setting No_TrainingNetworks to 10')
    No_TrainingNetworks = 10
  }
  
  errorMessage = packageMissingError("forecast", length(DataVec), SplitAt)
  if(errorMessage != FALSE){
    return(
      list(
        Model = NULL,
        Forecast = rep(NaN, length(DataVec)-SplitAt),
        TrainingData = DataVec,
        Info = errorMessage
      )
    )
  }
  
  errorMessage = packageMissingError("nnfor", length(DataVec), SplitAt)
  if(errorMessage != FALSE){
    return(
      list(
        Model = NULL,
        Forecast = rep(NaN, length(DataVec)-SplitAt),
        TrainingData = DataVec,
        Info = errorMessage
      )
    )
  }
  
  errorMessage = packageMissingError("DatabionicSwarm", length(DataVec), SplitAt)
  if(errorMessage != FALSE){
    return(
      list(
        Model = NULL,
        Forecast = rep(NaN, length(DataVec)-SplitAt),
        TrainingData = DataVec,
        Info = errorMessage
      )
    )
  }
  
  

  
  Time=as.Date(Time,tz=tz)
  if(Scaled){
    requireNamespace('DatabionicSwarm')
  if(!is.null(Predictors))
      PredictorTS=ts(data=DataVisualizations::RobustNormalization(Predictors),start=min(Time))
  else
    PredictorTS=NULL
  #A simple scaling procedure after Milligan and Cooper 1977
  #Save for later processing
  quants = quantile(DataVec, c(0.01, 0.5, 0.99), na.rm = F)
  minX = quants[1]
  maxX = quants[3]
  Denom = maxX - minX
  ResponseTS=ts(data=DataVisualizations::RobustNormalization(DataVec),start=min(Time))
  }else{
    if(!is.null(Predictors))
      PredictorTS=ts(data=Predictors,start=min(Time))
    else
      PredictorTS=NULL
    
    ResponseTS=ts(data=DataVec,start=min(Time))
    
   
  }
  #How many Weeks in Future should be predicted
  N=length(ResponseTS)
  if(!is.null(Predictors))
    if(nrow(PredictorTS)!=N) stop('Number of Cases of Predictors has to be equal to number of cases (length) of Datavectors.')
     
    
    TrainData=head(ResponseTS,SplitAt)
    TestData=tail(ResponseTS,N-SplitAt)
    ForecastTime=tail(Time,N-SplitAt)
  if(!is.null(Predictors)){
    TrainRegr=head(PredictorTS,SplitAt)
    TestRegr=tail(PredictorTS,N-SplitAt)
  }else{
    TrainRegr=NULL
    TestRegr=NULL
  }
  #print(str(TrainRegr))
  #Build model wit Trainings data
  if(PlotEvaluation){
    def.par <-
      par(no.readonly = TRUE) # save default, for resetting...
    m = graphics::layout(matrix(c(1, 2, 1, 2), 2, 2))
  }

  
  modelScaled=nnfor::elm(TrainData,outplot = PlotEvaluation,reps = No_TrainingNetworks,xreg = TrainRegr,hd=No_HiddenLayers,...)#,xreg.lags=c(1,1),direct = F)

  if(PlotEvaluation){
    title(xlab='Time in Unixtime units since 1970-01-01',ylab='Datavector',main=paste('Model on Training Data with',No_TrainingNetworks,'Reptitions'))#Underfitting?
    plot(modelScaled)
    par(def.par)
  }
  #Forecast with model on the part of Testdata of Predictors
  if(ForecastHorizon > N-SplitAt) {
    message("If predictors are used, forecast is only done on test data")
    predicted=forecast::forecast(modelScaled,h=N-SplitAt,xreg=PredictorTS)$mean
  } else {
    predicted=forecast::forecast(modelScaled,h=ForecastHorizon,xreg=PredictorTS)$mean
  }
  
  if(Scaled){ #Scale back
    predicted=predicted*Denom+minX
    TrainData=TrainData*Denom+minX
    TestData=TestData*Denom+minX
  }
  
  if(ForecastHorizon > N-SplitAt) {
    AccuracyTest=forecast::accuracy(predicted[1:length(TestData)],TestData)
  } else if(ForecastHorizon < N-SplitAt) {
    AccuracyTest=forecast::accuracy(predicted,TestData[1:ForecastHorizon])
  } else {
    AccuracyTest=forecast::accuracy(predicted,TestData)
  }
  
  if(PlotIt) {
    if(ForecastHorizon > N-SplitAt) {
      message("Plotting future evaluation/forecasts currently just implemented for the size of test data")
      predictedForPlot=predicted[1:length(TestData)]
    }
    get('plotEvaluationFilteredTS', envir = asNamespace('TSAT'), inherits = FALSE)(Time[(SplitAt+1):N],TestData,predicted[1:length(TestData)],FALSE)
  }
  
return(list(Model=modelScaled,Forecast=as.numeric(predicted),ForecastTime=ForecastTime,TrainingData=as.numeric(TrainData),TestData=as.numeric(TestData),Accuracy=AccuracyTest))
}