# DISC.R - Linear and quadratic discriminants for binary classification. # FIND LINEAR DISCRIMINANT. Takes as arguments a matrix of values for inputs in # a set of training cases (one case per row) and a vector of 0/1 class labels. # Returns a list in which element w is a vector pointing in the direction to # project inputs. Element w0 in the list is a scalar to add to the projected # value, chosen so that the classification error on the training cases is # minimized (when classifying to class 1 if the w %*% x + w0 > 0). lin.disc = function (X, t) { X0 = X[t==0,] X1 = X[t==1,] n0 = sum(t==0) n1 = sum(t==1) n = length(t) # Find means and sample covariance matrices for each class. m0 = apply(X0,2,mean) m1 = apply(X1,2,mean) v0 = cov(X0) v1 = cov(X1) # Find pooled covariance estimate. v = (n0/n)*v0 + (n1/n)*v1 # Find and return w and w0, which define the discriminant function. w = solve (v, m1-m0) d = as.vector (X %*% w) list (w=w, w0=find.w0(d,t)) } # FIND QUADRATIC DISCRIMINANT. Takes as arguments a matrix of values for # inputs in a set of training cases (one case per row) and a vector of 0/1 # class labels. Returns a list in which elements v (a matrix) and w (a vector) # define a quadratic function of an input vector, x, as x %*% v %*% x + w %*% x. # Element w0 in the list is a scalar to add to this value, chosen so that the # classification error on the training cases is minimized (when classifying to # class 1 if x %*% v %*% x + w %*%x > 0). quad.disc = function (X, t) { X0 = X[t==0,] X1 = X[t==1,] n0 = sum(t==0) n1 = sum(t==1) n = length(t) # Find means and sample covariance matrices for each class. m0 = apply(X0,2,mean) m1 = apply(X1,2,mean) v0 = cov(X0) v1 = cov(X1) # Find and return v, w and w0, which define the discriminant function. v = (solve(v0) - solve(v1)) / 2 w = solve(v1,m1) - solve(v0,m0) d = apply (X, 1, function (x) x %*% v %*% x + w %*% x) list (v=v, w=w, w0=find.w0(d,t)) } # FIND CONSTANT TERM FOR DISCRIMINANT THAT MINIMIZES CLASSIFICATION ERROR. # Given a vector of values for a discriminant function (without constant term) # for a set of training cases, and the vector of class labels for those cases, # this function returns the constant that should be added to the discriminant # so as to minimize the classification error rate on training cases. # # The same error rate will be produced by a range of values, from which the # midpoint is chosen. It's also possible that more than one range will produce # the same error rate, in which case one is picked arbitrarily. Values that # would classify all cases the same way are not considered. find.w0 = function (d, t) { # Reorder cases so that d will be in increasing order. o = order(d) d = d[o] t = t[o] # Choose amongst midpoints of adjacent values. Assume the midpoint of the # first range to start, then change if a better one is found. w0 = - (d[1] + d[2]) / 2 if (length(t)>1) { for (i in 2:length(t)) { w0.n = - (d[i-1] + d[i]) / 2 if (sum((d+w0.n>0)!=t) < sum((d+w0>0)!=t)) { w0 = w0.n } } } w0 }