/
oos_LINPROJ.R
executable file
·105 lines (99 loc) · 3.56 KB
/
oos_LINPROJ.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#' OOS : Linear Projection
#'
#' The simplest way of out-of-sample extension might be linear regression even though the original embedding
#' is not the linear type by solving
#' \deqn{\textrm{min}_{\beta} \|X_{old} \beta - Y_{old}\|_2^2} and use the estimate \eqn{\hat{beta}} to acquire
#' \deqn{Y_{new} = X_{new} \hat{\beta}}.
#'
#' @param Xold an \eqn{(n\times p)} matrix of data in original high-dimensional space.
#' @param Yold an \eqn{(n\times ndim)} matrix of data in reduced-dimensional space.
#' @param Xnew an \eqn{(m\times p)} matrix for out-of-sample extension.
#'
#' @return an \eqn{(m\times ndim)} matrix whose rows are embedded observations.
#'
#' @examples
#' \donttest{
#' ## generate sample data and separate them
#' data(iris, package="Rdimtools")
#' X = as.matrix(iris[,1:4])
#' lab = as.factor(as.vector(iris[,5]))
#' ids = sample(1:150, 30)
#'
#' Xold = X[setdiff(1:150,ids),] # 80% of data for training
#' Xnew = X[ids,] # 20% of data for testing
#'
#' ## run PCA for train data & use the info for prediction
#' training = do.pca(Xold,ndim=2)
#' Yold = training$Y
#' Ynew = Xnew%*%training$projection
#' Yplab = lab[ids]
#'
#' ## perform out-of-sample prediction
#' Yoos = oos.linproj(Xold, Yold, Xnew)
#'
#' ## visualize
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(1,2))
#' plot(Ynew, pch=19, col=Yplab, main="true prediction")
#' plot(Yoos, pch=19, col=Yplab, main="OOS prediction")
#' par(opar)
#' }
#'
#' @author Kisung You
#' @rdname oos_LINPROJ
#' @export
oos.linproj <- function(Xold, Yold, Xnew){
methodname = "linproj"
#------------------------------------------------------------------------
## PREPROCESSING
# 0. if inputs are vectors, make it a matrix.
if (is.vector(Xnew)){ Xnew = matrix(Xnew, nrow=1) }
if (is.vector(Xold)){ Xold = matrix(Xold, nrow=1) }
if (is.vector(Yold)){ Yold = matrix(Yold, nrow=1) }
# 1. parameters
n = nrow(Xold)
p = ncol(Xold)
ndim = ncol(Yold)
m = nrow(Xnew)
if (nrow(Yold)!=n){
stop(paste(" oos.",methodname," : 'Yold' should have same number of rows/observations as 'Xold'.",sep=""))
}
if (ndim >= p){
stop(paste("* oos.",methodname," : non-sense alert : reduced dimensionality should've been smaller than original data dimension."))
}
if (ncol(Xnew)!=p){
stop(paste("* oos.",methodname," : 'Xnew' should have same number of columns as 'Xold'.",sep=""))
}
# # 2. trfinfo
# if (!is.list(trfinfo)){
# stop(paste("* oos.",methodname," : 'trfinfo' should be provided a list.",sep=""))
# }
# if ((!('algtype'%in%names(trfinfo)))||(!('type'%in%names(trfinfo)))||(!('mean'%in%names(trfinfo)))||(!('multiplier'%in%names(trfinfo)))){
# stop(paste("* oos.",methodname," : 'trfinfo' is an invalid one. Use 'info' output you acquired from one of functions in the package."))
# }
# #------------------------------------------------------------------------
# # COMPUTE : PREPROCESSING
# # old ones
# X = aux.oospreprocess(Xold, trfinfo)
# Y = Yold
# # new ones
# XX = aux.oospreprocess(Xnew, trfinfo)
#
# #------------------------------------------------------------------------
# # COMPUTE : MAIN PART FOR LINEAR PROJECTION
# # 1. X\beta = Y
# betahat = aux.pinv(X)%*%Y
# # 2. use the projection matrix for this
# YY = XX%*%betahat
#
#
# #------------------------------------------------------------------------
# ## RETURN !
# result = list()
# result$Ynew = YY;
# return(result)
#------------------------------------------------------------------------
## COMPUTE
Ynew = oos_linproj(Xold, Yold, Xnew)
return(Ynew)
}