Demo

rm(list=ls())
setwd("F:/Working Section/Projects/FuncQuanReg/Codebase/R_code/QRFD_Rcode")
source('QuantileReg_FDA.R')

##############################################
# Model Settings
##############################################
seed = 790                          # Seed for random number generation
m <- 30                             # Number of grid points
n <- 1100                           # Sample size (1000 for training, 100 for testing)
sigma2 = 5^2                        # Variance of Y given X
eps.sigma2 = 0.25                   # Variance of noise for each trajectory
K = 4                               # Number of principle components
tau <- c(0.05,0.1,0.25,0.5,0.75,0.9,0.95)         # Quantile levels

ti <- seq(0,10,length=m)                  # Equally spaced design points on [0,10]         
t.seq <- t(matrix(rep(ti,n),m,n))

mu <- function(ti){ti + sin(ti)}                                # Mean function
mu.mat <- t(apply(t.seq,1,mu))

phi_1 = function(ti){cos(2*pi*ti/10)/sqrt(5)}                   # Eigenfunctiions
phi_2 = function(ti){sin(2*pi*ti/10)/sqrt(5)}
phi_3 = function(ti){cos(4*pi*ti/10)/sqrt(5)}
phi_4 = function(ti){sin(4*pi*ti/10)/sqrt(5)}
phi.mat = cbind(phi_1(ti), phi_2(ti), phi_3(ti), phi_4(ti))

lambda_1 = 4^2                                                  # Eigenvalues
lambda_2 = 3^2
lambda_3 = 2.75^2
lambda_4 = 2.25^2
Lambdasqrt.mat = diag(sqrt(c(lambda_1, lambda_2, lambda_3, lambda_4)))

##############################################
# Data Generation
##############################################
set.seed(seed)
xi.mat <- Lambdasqrt.mat %*% matrix(rnorm(4*n,mean=0,sd=1),nrow=4)
x.mat <- mu.mat + t(phi.mat %*% xi.mat) + 
  sqrt(eps.sigma2) * matrix(rnorm(n*m, mean=0, sd=1),nrow=n)    # Functional covariate
U <- runif(n,-16,16)                                            # Scalar covariate
Y <- rnorm(n,0,sqrt(sigma2)) + 2 * colSums(xi.mat) + 2 * U      # Response (normal)

# Split training/testing set
id.test <- sample(seq_len(n), size=100, replace=FALSE)
Y.test <- Y[id.test]; Y.training <- Y[!(seq_len(n) %in% id.test)]
U.test <- U[id.test]; U.training <- U[!(seq_len(n) %in% id.test)] 
X.test <- x.mat[id.test,]; X.training <- x.mat[!(seq_len(n) %in% id.test),]
scr <- cbind(xi.mat[,!(seq_len(n) %in% id.test)],xi.mat[,id.test])
data <- rbind(cbind(Y.training,U.training,X.training), cbind(Y.test,U.test,X.test))
id.test <- 1001:n

##############################################
# Joint Quantile Regression
##############################################
result <- QuantileReg_FDA(data,ti,id.test,tau,pve=0.99)

## Warning: package 'refund' was built under R version 3.2.4

##############################################
# Compare the predicted results with the truth
##############################################
predF <- result$predF.mono
predQ <- result$predQ
fixed.y <- result$fixed.y

mu <- 2 * colSums(scr) + 2 * data[,2]  
mu.test <- mu[id.test]
Fy <- t(sapply(mu.test, function(a) pnorm(fixed.y,a,5))) # Ture dis func for testing set
matplot(fixed.y,t(Fy),type="l",lwd=2,xlab="y",ylab="",main="True CDF")

matplot(fixed.y,t(predF),type="l",lwd=2,xlab="y",ylab="",main="Predicted CDF")

Qy <- matrix(NA,length(Y.test),length(tau))              # Ture quantile for testing set    
for (i in seq_len(length(Y.test)))
{
  Qy[i,] <- sapply(tau,function(a) ifelse(tail(Fy[i,],1) <= a, 
                                             tail(fixed.y,1),
                                             fixed.y[head(which(Fy[i,] > a),1)]))
}
plot(colMeans(Qy),colMeans(predQ),lwd=2,col="blue",xlim=c(-6,12),ylim=c(-6,12),
     xlab="True quantile",ylab="Predicted quantile",main="Q-Q plot")
abline(0,1,lwd=2,col="red")