############### Linear mFPCA functions ###############
#Use 'folds'-fold cross validation with 1-se rule to determine optimal number of principal components from 1 to "NumTestK" for local linear regression of y on x
K.Select.Lin=function(x,y,NumTestK=5,folds=10,PCA=F){
  n=nrow(x)
  ss=sample(rep(1:folds,diff(floor(n*seq(0,1,1/folds)))))
  Res=array(0,dim=c(ceiling(n/folds),folds,NumTestK))
  for(j in 1:folds){
    x1=x[which(ss!=j),]
    x2=x[which(ss==j),]
    y1=y[which(ss!=j)]
    y2=y[which(ss==j)]
    n1=nrow(x1)
    n2=nrow(x2)
    
    if(PCA==F){
      xmean=apply(x1,2,mean)
      x1center=x1-matrix(1,nrow=n1,ncol=1)%*%matrix(xmean,nrow=1)
      x2center=x2-matrix(1,nrow=n2,ncol=1)%*%matrix(xmean,nrow=1)
      
      #Eigen Decomposition
      eig=eigen(t(x1center)%*%x1center/n1) 
    }
    for(i in 1:NumTestK){
      if(PCA==F){
        #Eigenfunctions and Eigenvalues
        phi.hat=t(eig$vectors[,1:i])
        
        #Scores
        S1=x1center%*%t(phi.hat)
        S2=x2center%*%t(phi.hat)
        
        Yhat=try(S2%*%solve(t(S1)%*%S1)%*%t(S1)%*%y1,silent=T)
      }else{
        Yhat=try(x2%*%solve(t(x1)%*%x1)%*%t(x1)%*%y1,silent=T)
      }
      
      if(inherits(Yhat,"try-error")){
        Res[,j,i]=NA
      }else{
        Res[1:n2,j,i]=(y2-Yhat)^2/n2
      }
    }
  }
  CV=apply(Res,c(2,3),sum)
  CVest=apply(CV,2,mean)
  min=which.min(CVest)
  SE=apply(CV,2,sd)/sqrt(folds)
  InRange=CVest<(CVest[min]+SE[min])
  K=min(which(InRange))
  K
}

#Construct 1-alpha prediction interval for 'X.New' given (multiple) functional observations 'X' and response Y
#'X' is n by d matrix, with each row representing a single observation of (possibly multiple concatenated) functional predictors
#'X.New' is vector of length d
#'Y' is a vector of length n
#Number of (m)FPCA scores 'k' can be specified - otherwise cross-validation will select optimal values
#If PCA argument is True, function will treat 'X' and 'X.New' as functional principal component scores
#If contiguity condition is not met, algorithm will fail with the error message "X.New is outside the range of X"
lin.mFPCA=function(X,Y,X.New,k=NULL,PCA=F,alpha=.05){
  if(alpha <= 0 | alpha >= 1){
    stop("alpha must be within (0,1)")
  }
  if(!is.matrix(X)){
    stop("X must be an nxm matrix")
  }
  if(!is.matrix(X.New) & !is.vector(X.New)){
    stop("X.New must be a 1xm matrix or vector")
  }
  if(!is.matrix(Y) & !is.vector(Y)){
    stop("Y must be an nx1 matrix or vector")
  }
  n=nrow(X)
  
  require(Matrix)
  
  if(PCA){
    k=ncol(X)
  }
  
  if(is.null(k)){
    k=K.Select.Lin(X,Y)
  }
  
  if(n<2/alpha-1){
    stop("Sample size is too small for this alpha. Sample size must be at least 2/alpha - 1")
  }
  
  if(!PCA){
    #Center the n+1 observations
    X.Aug=rbind(X,X.New)
    X.mean=apply(X.Aug,2,mean)
    X.Aug.center=X.Aug-matrix(1,nrow=n+1,ncol=1)%*%matrix(X.mean,nrow=1)
    
    #Eigen Decomposition
    eig.n1=eigen(t(X.Aug.center)%*%X.Aug.center/(n+1)) 
    
    #Eigenfunctions and Eigenvalues
    phi.hat.n1=t(eig.n1$vectors[,1:min(n+1,k)])
    
    #Scores
    Scores.n1=X.Aug.center%*%t(phi.hat.n1)
    Scores.n=Scores.n1[1:n,]
  }else{
    Scores.n1=X.New
    Scores.n=X
  }
  
  #Fitting model first n obs.
  model.n=lm(Y~Scores.n)
  residuals.n=residuals(model.n)
  
  #Calculating Hat Matrix
  Score.Matrix.n1=cbind(rep(1,n+1),Scores.n1)
  Hat.Matrix=Score.Matrix.n1%*%solve(t(Score.Matrix.n1)%*%Score.Matrix.n1)%*%t(Score.Matrix.n1)
  
  #Test if X_n+1 is in range of X(n)
  if(sum((1-Hat.Matrix[n+1,n+1]+Hat.Matrix[1:n,n+1])<=0)>0){
    stop("X.New is outside the range of X")
  }
  
  if(sum((1-Hat.Matrix[n+1,n+1]-Hat.Matrix[1:n,n+1])<=0)>0){
    stop("X.New is outside the range of X")
  }
  
  #Conformity Scores
  a1=residuals.n/(1-Hat.Matrix[n+1,n+1]+Hat.Matrix[1:n,n+1])
  a2=-residuals.n/(1-Hat.Matrix[n+1,n+1]-Hat.Matrix[1:n,n+1])
  a=as.vector(sort(c(a1,a2),index.return=F))
  
  #Determine upper and lower bounds for residuals
  Lower=floor((n+1)*alpha)
  Upper=ceiling((2-alpha)*(n+1)-1)
  
  #Determine point estimate of y for the new observation using original n points
  y.n1.hat.n=as.numeric(model.n$coefficients[1]+Scores.n1[n+1,]%*%model.n$coefficients[-1])
  
  #Determine Conformal Prediction Interval
  Y.low=y.n1.hat.n+a[Lower]
  Y.high=y.n1.hat.n+a[Upper]
  
  list("LB"=Y.low,"UB"=Y.high)
}

############### Local Linear mFPCA functions ###############
# Obtain first k (m)FPCA scores of x
get.PC.Scores=function(x,xnew,k=NULL){
  n=nrow(x)
  X.Aug=rbind(x,xnew)
  X.mean=apply(X.Aug,2,mean)
  X.Aug.center=X.Aug-matrix(1,nrow=nrow(X.Aug),ncol=1)%*%matrix(X.mean,nrow=1)
  
  #Eigen Decomposition
  eig.n1=eigen(t(X.Aug.center)%*%X.Aug.center/(n+1)) 
  
  #Eigenfunctions and Eigenvalues
  phi.hat.n1=t(eig.n1$vectors[,1:min(n+1,k)])
  
  #Scores
  Scores.n1=X.Aug.center%*%t(phi.hat.n1)
  Scores.n1
}

#Construct X matrix based on FPCA scores at point 'xnew'
create.X.matrix=function(x,xnew){
  if(is.vector(xnew)){
    if(ncol(x)==length(xnew)){
      n=nrow(x)
      p=ncol(x)
      X1=matrix(1,n,p)
      for(i in 1:n){
        X1[i,]=x[i,]-xnew
      }
      X=cbind(rep(1,n),X1)
      X
    }else{
      stop("x and xnew must have same number of columns (create.X.matrix)")
    }
  }else{
    stop("xnew must be vector (create.X.matrix)")
  }
}

# Construct weight matrix based on FPCA scores at point 'xnew', bandwidth h, and specified kernel
create.W.matrix=function(x,xnew,h=NULL,Kernel="Asym.Norm"){
  if(is.vector(xnew)){
    if(ncol(x)==length(xnew)){
      dists=rep(0,nrow(x))
      for(i in 1:nrow(x)){
        dists[i]=sqrt(sum(t(t(x[i,]-xnew)/h)^2))
      }
      if(Kernel=="Asym.Norm"){
        weights=2*dnorm(dists)
        W=diag(weights)
        W
      }else if(Kernel=="Epanechnikov"){
        weights=1.5*(1-(dists)^2)
        weights=ifelse(weights<0,0,weights)
        W=diag(weights)
        W
      }else{
        stop("Unrecognized kernel (create.W.matrix)")
      }
    }else{
      stop("x and xnew must have same number of columns (create.W.matrix)")
    }
  }else{
    stop("xnew must be a vector (create.W.matrix)")
  }
}

#Estimate the beta for the local linear regression of y on x
Create.Beta = function(x,xnew,y,...){
  if(is.vector(xnew)){
    if(ncol(x)==length(xnew)){
      if(nrow(x)==(length(y)+1)){
        n=nrow(x)
        X=create.X.matrix(x,xnew)
        Xn=X[1:(n-1),]
        W=create.W.matrix(x,xnew,...)
        Wn=W[1:(n-1),1:(n-1)]
        Beta=solve(t(Xn)%*%Wn%*%Xn)%*%t(Xn)%*%Wn%*%y
        Beta
      }else if(nrow(x)==(length(y))){
        n=nrow(x)
        X=create.X.matrix(x,xnew)
        W=create.W.matrix(x,xnew,...)
        Beta=solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%y
        Beta
      }
      else{
        stop("y must be of length n (Create.Beta)")
      }
    }else{
      stop("x and xnew must have same number of columns (Create.Beta)")
    }
  }else{
    stop("xnew must be a vector (Create.Beta)")
  }
}

#Calculate the A_i from the algorithm
Create.A = function(x,xnew,y,...){
  if(is.vector(xnew)){
    if(ncol(x)==length(xnew)){
      n=nrow(x)
      p=ncol(x)
      e=c(1,rep(0,p))
      X=create.X.matrix(x,xnew)
      Xn=X[1:n-1,]
      Xn1=X[n,]
      W=create.W.matrix(x,xnew,...)
      Wn=W[1:(n-1),1:(n-1)]
      Wn1=W[n,n]
      Beta=Create.Beta(x,xnew,y,...)
      I=diag(1,p+1)
      A=t(e)%*%(I-solve(t(Xn)%*%Wn%*%Xn)%*%Xn1%*%Wn1%*%t(Xn1)/(1+as.numeric(t(Xn1)%*%solve(t(Xn)%*%Wn%*%Xn)%*%Xn1*Wn1)))%*%Beta
      as.numeric(A)
    }else{
      stop("x and xnew must have same number of columns (create.A)")
    }
  }else{
    stop("xnew must be a vector (create.A)")
  }
}

#Calculate the B_i from the algorithm
Create.B = function(x,xnew,...){
  if(ncol(x)==length(xnew)){
    if(is.vector(xnew)){
      n=nrow(x)
      p=ncol(x)
      e=c(1,rep(0,p))
      X=create.X.matrix(x,xnew)
      Xn=X[1:n-1,]
      Xn1=X[n,]
      W=create.W.matrix(x,xnew,...)
      Wn=W[1:(n-1),1:(n-1)]
      Wn1=W[n,n]
      B=t(e)%*%solve(t(Xn)%*%Wn%*%Xn)%*%Xn1%*%Wn1/(1+as.numeric(t(Xn1)%*%solve(t(Xn)%*%Wn%*%Xn)%*%Xn1*Wn1))
      as.numeric(B)
    }else{
      stop("xnew must be a vector (Create.B)")
    }
  }else{
    stop("x and xnew must have same number of columns (Create.B)")
  }
}

#Use 'folds'-fold cross validation with 1-se rule to determine optimal number of principal components from 1 to "NumTestK" for local linear regression of y on x
#One can specify bandwith 'h' to use
#Non-functional covariates can be included by including them in the "Covar" argument
k.Select.ll=function(x,y,NumTestK=5,folds=10,h=NULL,Covar=NULL,PCA=F,...){
  nullh=is.null(h)
  n=nrow(x)
  ss=sample(rep(1:folds,diff(floor(n*seq(0,1,1/folds)))))
  Res=array(0,dim=c(ceiling(n/folds),folds,NumTestK))
  for(j in 1:folds){
    x1=x[which(ss!=j),]
    x2=x[which(ss==j),]
    y1=y[which(ss!=j)]
    y2=y[which(ss==j)]
    n1=nrow(x1)
    n2=nrow(x2)
    xmean=apply(x1,2,mean)
    x1center=x1-matrix(1,nrow=n1,ncol=1)%*%matrix(xmean,nrow=1)
    x2center=x2-matrix(1,nrow=n2,ncol=1)%*%matrix(xmean,nrow=1)
    if(!is.null(Covar)){
      Covar1=Covar[which(ss!=j),]
      Covar2=Covar[which(ss==j),]
    }
    
    #Eigen Decomposition
    eig=eigen(t(x1center)%*%x1center/n1) 
    for(i in 1:NumTestK){
      #Eigenfunctions and Eigenvalues
      phi.hat=t(eig$vectors[,1:i])
      
      #Scores
      if(!PCA){
        Scores1=x1center%*%t(phi.hat)
        Scores2=x2center%*%t(phi.hat)
      }else{
        Scores1=x1
        Scores2=x2
      }
      
      if(!is.null(Covar)){
        Scores1=cbind(Scores1,Covar1)
        Scores2=cbind(Scores2,Covar2)
      }
      
      if(nullh){
        sd=apply(Scores1,2,sd)
        h=5*sd*n1^(-.2)
      }
      for(k in 1:n2){
        X=create.X.matrix(Scores1,Scores2[k,])
        W=create.W.matrix(Scores1,Scores2[k,],h)
        
        Yhat=try(solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%y1,silent=T)
        if(inherits(Yhat,"try-error")){
          Res[k,j,i]=NA
        }else{
          Res[k,j,i]=(y2[k]-Yhat[1])^2/n2
        }
      }
    }
  }
  CV=apply(Res,c(2,3),sum)
  CVest=apply(CV,2,mean)
  min=which.min(CVest)
  SE=apply(CV,2,sd)/sqrt(folds)
  InRange=CVest<(CVest[min]+SE[min])
  K=min(which(InRange))
  K
}

#Use 'folds'-fold cross validation with 1-se rule to determine optimal number of principal components from 1 to "NumTestH" for local linear regression of y on x
#Non-functional covariates can be included by including them in the "Covar" argument
h.Select.ll=function(x,y,NumTestH=20,folds=10,K,Covar=NULL,PCA=F,...){
  n=nrow(x)
  HTest=matrix(0,nrow=NumTestH,ncol=K)
  if(!PCA){
    xmean=apply(x,2,mean)
    xcenter=x-matrix(1,nrow=n,ncol=1)%*%matrix(xmean,nrow=1)
    eig=eigen(t(xcenter)%*%xcenter/n)
    phi.hat=t(eig$vectors[,1:K])
    Scores=xcenter%*%t(phi.hat)
  } else{
    Scores=x
  }
  if(!is.null(Covar)){
    Scores=cbind(Scores,Covar)
  }
  
  for(i in 1:K){
    HTest[,i]=exp(seq(log(max(abs(dist(Scores[,i])/40))),log(max(abs(dist(Scores[,i])/2))),length.out=NumTestH))
  }
  ss=sample(rep(1:folds,diff(floor(n*seq(0,1,1/folds)))))
  Res=array(0,dim=c(ceiling(n/folds),folds,NumTestH))
  for(j in 1:folds){
    x1=x[which(ss!=j),]
    x2=x[which(ss==j),]
    y1=y[which(ss!=j)]
    y2=y[which(ss==j)]
    n1=nrow(x1)
    n2=nrow(x2)
    if(!is.null(Covar)){
      Covar1=Covar[which(ss!=j),]
      Covar2=Covar[which(ss==j),]
    }
    if(!PCA){
      xmean=apply(x1,2,mean)
      x1center=x1-matrix(1,nrow=n1,ncol=1)%*%matrix(xmean,nrow=1)
      x2center=x2-matrix(1,nrow=n2,ncol=1)%*%matrix(xmean,nrow=1)
      
      
      #Eigen Decomposition
      eig=eigen(t(x1center)%*%x1center/(n1)) 
      #Eigenfunctions and Eigenvalues
      phi.hat=t(eig$vectors[,1:K])
      
      #Scores
      Scores1=x1center%*%t(phi.hat)
      Scores2=x2center%*%t(phi.hat)
    } else{
      Scores1=x1
      Scores2=x2
    }
    if(!is.null(Covar)){
      Scores1=cbind(Scores1,Covar1)
      Scores2=cbind(Scores2,Covar2)
    }
    
    for(i in 1:NumTestH){
      for(k in 1:n2){
        X=create.X.matrix(Scores1,Scores2[k,])
        W=create.W.matrix(Scores1,Scores2[k,],HTest[i,])
        Yhat=try(solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%y1,silent=T)
        if(inherits(Yhat,"try-error")){
          Res[k,j,i]=NA
        }else{
          Res[k,j,i]=(y2[k]-Yhat[1])^2/n2
        }
      }
    }
  }
  CV=apply(Res,c(2,3),sum)
  CVest=apply(CV,2,mean)
  min=which.min(CVest)
  SE=apply(CV,2,sd)/sqrt(folds)
  InRange=CVest<(CVest[min]+SE[min])
  H=HTest[max(which(InRange)),]
  H
}

#Construct 1-alpha prediction interval for 'xnew' given (multiple) functional observations 'x' and response y
#'x' is n by d matrix, with each row representing a single observation of (possibly multiple concatenated) functional predictors
#'xnew' is vector of length d
#'y' is a vector of length n
#bandwidth 'h' and number of (m)FPCA scores 'k' can be specified - otherwise cross-validation will select optimal values
#If PCA argument is True, function will treat 'x' and 'xnew' as functional principal component scores
#If contiguity condition is not met, algorithm will fail with the error message "xnew must be within range of x"
ll.mFPCA=function(x,xnew,y,alpha=0.05,Kernel="Asym.Norm",h=NULL,k=NULL,PCA=F,...){
  if(alpha>0 & alpha<1){
    if(!PCA){
      if(is.null(k)){
        k=k.Select.ll(x,y,...)
      }
      Scores=get.PC.Scores(x,xnew,k,...)
    } else{
      if(is.null(k)){
        k=ncol(x)
      }
      Scores=rbind(x,xnew)
    }
    
    if(is.null(h)){
      h=h.Select.ll(x,y,K=k,PCA=PCA,...)
    }
    
    n=nrow(x)
    cond1=rep(0,n)
    cond2=rep(0,n)
    Bn1=Create.B(Scores,Scores[n+1,],h,Kernel,...)
    An1=Create.A(Scores,Scores[n+1,],y,h,Kernel,...)
    A=rep(0,n)
    B=rep(0,n)
    a=rep(0,n)
    b=rep(0,n)
    c=rep(0,n)
    d=rep(0,n)
    for(i in 1:n){
      B[i]=Create.B(Scores,Scores[i,],h,Kernel,...)
      A[i]=Create.A(Scores,Scores[i,],y,h,Kernel,...)
      cond1[i]=1-Bn1+B[i]
      cond2[i]=1-Bn1-B[i]
      a[i]=(y[i]-A[i]+An1)/(1-Bn1+B[i])
      b[i]=(-y[i]+A[i]+An1)/(1-Bn1-B[i])
      c[i]=min(a[i],b[i])
      d[i]=max(a[i],b[i])
    }
    if(sum(cond1<0)>0 | sum(cond2<0)>0){
      list("LB"=NA,"UB"=NA,"Bandwidth"=h)
      #stop("xnew must be within range of x")
    }else{
      c=sort(c)
      d=sort(d)
      Lower=floor((n+1)*alpha)
      Upper=ceiling((1-alpha)*(n+1))
      
      Y.low=c[Lower]
      Y.high=d[Upper]
      
      list("LB"=Y.low,"UB"=Y.high,"Bandwidth"=h)
    }
  }else{
    stop("alpha must be in (0,1)")
  }
}

############### Signature multiple Partial Linear functions ###############
h.SelectPL=function(Y,Z,NumTestH=20,folds=10){
  Y=as.matrix(Y)
  Z=as.matrix(Z)
  n=nrow(Y)
  K=ncol(Z)
  Max=rep(0,K)
  HTest=matrix(0,nrow=NumTestH,ncol=K)
  for(i in 1:K){
    Max[i]=max(abs(dist(Z[,i])))
    HTest[,i]=exp(seq(log(Max[i]/40),log(Max[i]/2),length.out=NumTestH))
  }
  ss=sample(rep(1:folds,diff(floor(n*seq(0,1,1/folds)))))
  Res=array(0,dim=c(ceiling(n/folds),folds,NumTestH))
  for(j in 1:folds){
    Y1=as.matrix(Y[which(ss!=j),])
    Y2=as.matrix(Y[which(ss==j),])
    Z1=Z[which(ss!=j),]
    Z2=Z[which(ss==j),]
    n1=nrow(Y1)
    n2=nrow(Y2)
    for(i in 1:NumTestH){
      Wm=Create.Wm(Z1,HTest[i,])
      mhats=Wm%*%Y1
      newmhats=matrix(0,nrow=n2,ncol=K)
      for(k in 1:K){
        mhat=mhats[((k-1)*n1+1):(k*n1),]
        for(l in 1:n2){
          if(all(Z1[,k]>Z2[l,k])){
            ClosestAbove=min(which(Z1[,k]==min(Z1[Z1[,k]>=Z2[l,k],k])))
            newmhats[l,k]=mhat[ClosestAbove]
          }else if(all(Z1[,k]<Z2[l,k])){
            ClosestBelow=min(which(Z1[,k]==max(Z1[Z1[,k]<=Z2[l,k],k])))
            newmhats[l,k]=mhat[ClosestBelow]
          }else{
            ClosestBelow=min(which(Z1[,k]==max(Z1[Z1[,k]<=Z2[l,k],k])))
            ClosestAbove=min(which(Z1[,k]==min(Z1[Z1[,k]>=Z2[l,k],k])))
            DistBelow=(Z2[l,k]-Z1[ClosestBelow,k])
            DistAbove=(Z1[ClosestAbove,k]-Z2[l,k])
            TotalDist=DistBelow+DistAbove
            WeightBelow=DistAbove/TotalDist
            WeightAbove=DistBelow/TotalDist
            newmhats[l,k]=WeightBelow*mhat[ClosestBelow]+WeightAbove*mhat[ClosestAbove]
          }
        }
      }
      Yhat=rowSums(newmhats)
      Res[,j,i]=(Y2-Yhat)^2/n2
    }
  }
  CV=apply(Res,c(2,3),sum)
  CVest=apply(CV,2,mean)
  min=which.min(CVest)
  SE=apply(CV,2,sd)/sqrt(folds)
  InRange=CVest<(CVest[min]+SE[min])
  H=HTest[max(which(InRange)),]
  H
}

Create.Wm=function(Z,h){
  Z=as.matrix(Z)
  n=nrow(Z)
  d=ncol(Z)
  S=matrix(0,nrow=n*d,ncol=n)
  for(i in 1:n){
    for(j in 1:d){
      dist=(Z[,j]-Z[i,j])/h[j]
      s=dnorm(dist)/sum(dnorm(dist))
      S[(j-1)*n+i,]=s
    }
  }
  cent=diag(rep(1,n))-matrix(1/n,nrow=n,ncol=n)
  for(j in 1:d){
    S[((j-1)*n+1):(j*n),]=cent%*%S[((j-1)*n+1):(j*n),]
  }
  M=matrix(rep(S,d),nrow=d*n)
  for(j in 1:d){
    M[((j-1)*n+1):(j*n),((j-1)*n+1):(j*n)]=diag(n)
  }
  Minv=solve(M)
  W=Minv%*%S+matrix(1/(d*n),nrow=n*d,ncol=n)
  W
}

mPL.Sig=function(X,Z,XNew,ZNew,Y,alpha=0.05,h=NULL){
  n=nrow(X)
  p=ncol(X)
  d=ncol(Z)
  Z=as.matrix(Z)
  XAug=as.matrix(rbind(X,XNew))
  ZAug=as.matrix(rbind(Z,ZNew))
  if(is.null(h)){
    h=rep(0,d)
    for(j in 1:d){
      h[j]=max(as.matrix(dist(Z[,j])))/10
    }
  }
  W=Create.Wm(ZAug,h)
  Wx=W%*%XAug
  Wy=W[,1:n]%*%Y
  gx=matrix(0,nrow=n+1,ncol=p)
  gyn=rep(0,n+1)
  gyn1=rep(0,n+1)
  for(j in 1:d){
    gx=gx+Wx[((j-1)*(n+1)+1):(j*(n+1)),]
    gyn=gyn+Wy[((j-1)*(n+1)+1):(j*(n+1)),]
    gyn1=gyn1+W[((j-1)*(n+1)+1):(j*(n+1)),n+1]
  }
  Rx=XAug-gx
  RxRx=solve(t(Rx)%*%Rx)
  RxY=c(Y)%*%Rx[1:n,]
  RxKYn=gyn%*%Rx
  RxKYn1=gyn1%*%Rx
  A=Rx%*%RxRx%*%t(RxKYn-RxY)-gyn
  B=Rx%*%RxRx%*%t(RxKYn1-Rx[n+1,])-gyn1
  
  if(sum(abs(B[1:n])>1+B[n+1])>0){
    list("LB"=NA,"UB"=NA)
  }else{
    
    Lim1=(Y+A[1:n]-A[n+1])/(1+B[n+1]-B[1:n])
    Lim2=(-Y-A[1:n]-A[n+1])/(1+B[n+1]+B[1:n])
    Min=pmin(Lim1,Lim2)
    Max=pmax(Lim1,Lim2)
    Min=sort(Min)
    Max=sort(Max)
    Lower=floor((n+1)*alpha)
    Upper=ceiling((1-alpha)*(n+1))
    
    LB=Min[Lower]
    UB=Max[Upper]
    
    list("LB"=LB,"UB"=UB)
  }
}
