forked from jkrijthe/RSSL
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathNormalBasedClassifier.R
More file actions
165 lines (134 loc) · 7.51 KB
/
Copy pathNormalBasedClassifier.R
File metadata and controls
165 lines (134 loc) · 7.51 KB
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
#' @include Classifier.R
setClass("NormalBasedClassifier",
representation(means="matrix",sigma="list",prior="matrix"),
prototype(name="NormalBased Classifier"),
contains="Classifier")
#' @rdname rssl-formatting
#' @aliases show,NormalBasedClassifier-method
setMethod("show", signature(object="NormalBasedClassifier"), function(object) {
cat(object@name,'\n')
cat("Means:\n")
print(object@means)
cat("\nSigma:\n")
print(object@sigma)
cat("\nScaling:\n")
print(object@scaling)
cat("Classes:")
print(object@classnames)
})
#' @rdname rssl-predict
#' @aliases predict,NormalBasedClassifier-method
setMethod("predict", signature(object="NormalBasedClassifier"), function(object,newdata) {
ModelVariables<-PreProcessingPredict(modelform=object@modelform,
newdata=newdata,
y=NULL,
scaling=object@scaling,
intercept=FALSE)
X <- ModelVariables$X
M <- object@means
G <- matrix(NA,nrow(X),length(object@sigma))
for (c in 1:length(object@sigma)) {
S <- object@sigma[[c]]
C <- -ncol(X) * log(2 * pi) * 0.5 - 0.5 * log(det(S))
G[,c] <- C + c(log(object@prior[c,,drop=FALSE])) - rowSums( ((X-matrix(1,nrow(X),1) %*% M[c,,drop=FALSE]) %*% solve(S)) * (X-matrix(1,nrow(X),1) %*% M[c,,drop=FALSE]))/2
}
factor(as.numeric(which_rowMax(G)), labels=object@classnames,levels=1:length(object@classnames))
})
#' Return the negative log likelihood on the given dataset
#' @rdname loss-methods
#' @aliases loss,NormalBasedClassifier-method
setMethod("loss", signature(object="NormalBasedClassifier"), function(object, newdata, y=NULL) {
ModelVariables<-PreProcessingPredict(object@modelform,newdata,y=y,object@scaling,intercept=FALSE,classnames=object@classnames)
X<-ModelVariables$X
y<-ModelVariables$y
if (is.null(y)) { stop("No labels supplied.")}
Y <- model.matrix(~as.factor(y)-1)
k<-ncol(X) # Number of features
m<-object@means
ll<-0
losses <- rep(NA,nrow(X))
for (c in 1:nrow(m)) {
sigma<-object@sigma[[c]]
Xc <- X[as.logical(Y[,c]), ,drop=FALSE] #Select all objects in class c
ll <- ll + nrow(Xc) * ( log(object@prior[c,])-(k/2)*log(2*pi)-(1/2)*log(det(sigma)) ) #Add the constant part for each row in Xc
ll <- ll + -(1/2)*sum((Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma) * (Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])) #Add the dynamic contribution
# Note: much faster than: sum(-(1/2)*diag((Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma)%*%t(Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])))
losses[as.logical(Y[,c])] <- c( log(object@prior[c,])-(k/2)*log(2*pi)-(1/2)*log(det(sigma)) ) + -(1/2)*rowSums((Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma) * (Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE]))
}
return(-losses)
})
#' @param newdata design matrix
#' @param Y class responsibility matrix
#' @rdname losspart-methods
#' @aliases losspart,NormalBasedClassifier-method
setMethod("losspart", signature(object="NormalBasedClassifier"), function(object, newdata, Y) {
ModelVariables<-PreProcessingPredict(object@modelform,newdata,y=NULL,object@scaling,intercept=FALSE)
X<-ModelVariables$X
k<-ncol(X) # Number of features
m<-object@means # Parameters of the NM classifier
ll<-0
for (c in 1:nrow(m)) {
sigma<-object@sigma[[c]]
Xc<-X
ll<-ll + sum(Y[,c]) * ( log(object@prior[c,])-(k/2)*log(2*pi)-(1/2)*log(det(sigma)) ) #Add the constant part for each row in Xc
ll<-ll + sum(t(Y[,c,drop=FALSE]) %*% (-(1/2)*(Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma) * (Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE]))) #Add the dynamic contribution
# Note: much faster than: sum(-(1/2)*diag((Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma)%*%t(Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])))
}
return(-ll)
})
#' @param newdata Design matrix of labeled objects
#' @param Y label matrix of labeled objects
#' @param X_u Design matrix of unlabeled objects
#' @param Y_u label matrix of unlabeled objects
#' @rdname losslogsum-methods
#' @aliases losslogsum,NormalBasedClassifier-method
setMethod("losslogsum", signature(object="NormalBasedClassifier"), function(object, newdata, Y, X_u, Y_u) {
ModelVariables<-PreProcessingPredict(object@modelform,newdata,y=NULL,object@scaling,intercept=FALSE)
X<-ModelVariables$X
k<-ncol(X) # Number of features
m<-object@means # Parameters of the NM classifier
ll<-0
sigma<-object@sigma[[1]]
for (c in 1:nrow(m)) {
sigma<-object@sigma[[c]]
Xc<-X
ll<-ll + sum(Y[,c]) * ( log(object@prior[c,])-(k/2)*log(2*pi)-(1/2)*log(det(sigma)) ) #Add the constant part for each row in Xc
ll<-ll + sum(t(Y[,c,drop=FALSE]) %*% (-(1/2)*(Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma) * (Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE]))) #Add the dynamic contribution
# Note: much faster than: sum(-(1/2)*diag((Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma)%*%t(Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])))
}
Xc<-X_u
for (i in 1:nrow(Y_u)) {
ll <- ll + log(exp(( log(object@prior[1,])-(k/2)*log(2*pi)-(1/2)*log(det(sigma)) ) + ( sum(-(1/2)*(Xc[i,]-matrix(1,1,1) %*% m[1, ,drop=FALSE])%*%solve(sigma) * (Xc[i,]-matrix(1,1,1) %*% m[1, ,drop=FALSE])))) + exp(( log(object@prior[2,])-(k/2)*log(2*pi)-(1/2)*log(det(sigma)) ) + ( sum(-(1/2)*(Xc[i,]-matrix(1,1,1) %*% m[2, ,drop=FALSE])%*%solve(sigma) * (Xc[i,]-matrix(1,1,1) %*% m[2, ,drop=FALSE]))))) #Add the constant part for each row in Xc
}
# Note: much faster than: sum(-(1/2)*diag((Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])%*%solve(sigma)%*%t(Xc-matrix(1,nrow(Xc),1) %*% m[c, ,drop=FALSE])))
return(-ll)
})
#' @param newdata matrix of dataframe of objects to be classified
#' @rdname posterior-methods
#' @aliases posterior,NormalBasedClassifier-method
setMethod("posterior", signature(object="NormalBasedClassifier"), function(object,newdata) {
ModelVariables<-PreProcessingPredict(modelform=object@modelform,
newdata=newdata,
y=NULL,
scaling=object@scaling,
intercept=FALSE)
X<-ModelVariables$X
M<-object@means
G <- matrix(NA,nrow(X),length(object@sigma))
for (c in 1:length(object@sigma)) {
S <- object@sigma[[c]] # Covariance matrix
k <- ncol(S) # Dimensionality
G[,c] <- c(log(object@prior[c,,drop=FALSE]) - (k/2) * log(2*pi) - 0.5*log(det(S))) - (1/2)*rowSums( ((X-matrix(1,nrow(X),1) %*% M[c,,drop=FALSE]) %*% solve(S)) * (X-matrix(1,nrow(X),1) %*% M[c,,drop=FALSE]))
}
posteriors <- G - (log(rowSums(exp(G-apply(G,1,max))))+apply(G,1,max)) # More stable way of doing logsumexp
posteriors <- exp(posteriors)
colnames(posteriors) <- object@classnames
return(posteriors)
})
#' @rdname line_coefficients-methods
#' @aliases line_coefficients,LeastSquaresClassifier-method
setMethod("line_coefficients",signature(object="NormalBasedClassifier"), function(object) {
w <- -(object@means[2,, drop=FALSE]-object@means[1,, drop=FALSE]) %*% solve(object@sigma[[1]])
w0 <-(log(object@prior[1])-log(object@prior[2])- 0.5*object@means[1,, drop=FALSE] %*% solve(object@sigma[[1]]) %*% t(object@means[1,, drop=FALSE]) + 0.5*object@means[2,, drop=FALSE] %*% solve(object@sigma[[1]]) %*% t(object@means[2,, drop=FALSE]))
return(coefficients_after_scaling(w0=w0,w=w,scaling=object@scaling))
})