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(xclass) = 1/sqrt(2pi^d * Sigma ) * exp( (xmu)Sigma^1(xmu))
>
> 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 ddimensional, 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.
Apparently classify takes the middle road.
Hope this helps.
Greg
