marquito wrote: > Hi! > > In my thesis I'm using to classify data the Matlab 'classify' > function with linear discrimination. To see what the math behind it > is, I looked into the code and found this: > > ===================================================== > % Pooled estimate of covariance > [Q,R] = qr(training - gmeans(gindex,:), 0); > R = R / sqrt(n - ngroups); % SigmaHat = R'*R > s = svd(R); > if any(s <= eps^(3/4)*max(s)) > error('The pooled covariance matrix of TRAINING > must be positive definite.'); > end > > % MVN relative log posterior density, by group, for > each sample > for k = 1:ngroups > A = (sample - repmat(gmeans(k,:), mm, 1)) / R; > D(:,k) = log(prior(k)) - .5*sum(A .* A, 2); > end > ====================================================== > > I dont know exactly, was is going on there. I expected to see > something like a multivariate Gauss distribution, like: > > p(x|class) = 1/sqrt(2pi^d * |Sigma| ) * exp( (x-mu)Sigma^-1(x-mu)) > > or something similar to this. Could somebody verify this or explain > what kind of magic the programmer used (why qr decomposition?). > I'm also very interested how 'D' is calculated since I use this value > to show the distances of a sample to the different classes. Shouldn't > this be the probability density function of x for the different > classes?
Just because one postulates that x is d-dimensional, doesn't make it so. Accordingly, sometimes two or more inputs are highly collinear and the data really lies in a lower dimensional subspace. Sigma is then ill conditioned (i.e., nearly singular) and, at best, matrix inversion results in garbage. However, you don't have to cry "error" and give up.
There are at least two fixes. First, reduce the dimensionality of the data by projecting it into the dominant subspace. Second, use matrix pseudoinversion. I successfully used the latter for more than 20 years.
QR decomposition yields more accurate results than ordinary inversion as long as the matrix is not too badly conditioned. The conditioning of sigma is readily determined by using rank(sigma) and cond(sigma). If it is well conditioned, the slash solution sigma\eye(d) will work. If it is not well conditioned but not quite singular, use QR. If it is ill conditioned or even singular, use pseudoinversion.