R/pcr.r

In QWDAP: Quantum Walk-Based Data Analysis and Prediction

#' @title Principle Component Regression
#' @description Principle component regression. This is a linear regression method used to
#' establish the linear relationship between the original time series and the modes
#' generated by quantum walks.
#' @usage qwdap.pcr(in_data, data_range, plotting)
#' @param in_data a 'QWMS' object, which includes the target series and the
#' selected modes which can be obtained from modes selection.
#' @param data_range the range of the train samples.
#' @param plotting whether to plot.
#'
#' @return a 'QWMODEL' object which includes the information of regression analysis.
#' @import pls
#' @importFrom stats as.formula
#' @importFrom graphics lines legend
#' @export qwdap.pcr
#'
#' @examples
#' data("traffic.n1")
#' res.pcr <- qwdap.pcr(traffic.n1,c(1,500), FALSE)
#' 
qwdap.pcr<-function(in_data, data_range, plotting = FALSE){
  # library(pls)
  if(class(in_data)!='QWMS'){
    stop("The 'in_data' is not a 'QWMS' object.")
  }
  if(!is.vector(data_range)||!is.numeric(data_range)||length(data_range)<2){
    stop("The parameter 'data_range' is error.")
  }
  # pre combine
  co_data = cbind(in_data$real, in_data$ctqw)
  co_data <- subset(co_data, select = c(colnames(in_data$real), in_data$variate))
  # data_y=in_data$real[data_range[1]:data_range[2],]
  # data_x=in_data$ctqw[data_range[1]:data_range[2],]
  #colnames(data_y)<-c("y")
  # lm_1 = colnames(data_y)
  # lm_2 = colnames(data_x)
  my_lm = paste(colnames(co_data)[1],"~",paste(colnames(co_data)[-1],collapse = " + "))
  # my_lm = lm(as.formula(my_lm),data = cbind(data_y,data_x))
  res<-pcr(as.formula(my_lm),data = as.data.frame(co_data[data_range[1]:data_range[2],]),validation = "CV")
  if(plotting){
    data_plot = as.data.frame(res$fitted.values)
    # par(mfrow=c(ceiling(ncol(data_x)/2.0),if(ncol(data_x)==1) 1 else 2))
    # tmp_data = cbind(data_y,data_plot)
    for(i in (1:ncol(data_plot))){
      # if(i==1){
      #   plot(x=c(1:nrow(tmp_data)),y=tmp_data[,1],type = "l",xlab="index",ylab="value",lwd=2)
      # }
      # lines(x=c(1:nrow(tmp_data)),y=tmp_data[,i],type = "l",col=i,lwd=2)
      tmp_data = cbind(in_data$real[data_range[1]:data_range[2],],data_plot[,i])
      colnames(tmp_data)=c("Actual series","Fitted series")
      plot(x=c(1:nrow(tmp_data)),y=tmp_data[,1],type = "l",xlab="index",ylab="value",
           ylim=c(min(tmp_data)-2,max(tmp_data)+2),lwd=1)
      lines(x=c(1:nrow(tmp_data)),y=tmp_data[,2],type = "l",col=2,lwd=1)
      legend("topleft", colnames(tmp_data), col = c(1,2),
             lwd = c(1), bg = "grey95", box.col = NA,
             cex = 0.8, inset = c(0.02, 0.03), ncol = 1)
    }
  }
  res<-list(real = in_data$real, ctqw = co_data[,-1], index = in_data$index,
            method = "PCR",model=res)
  res<-structure(res,class="QWMODEL")
  return(res)
}