help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: how can I do Principal Components Analysis with octave?


From: Stefan van der Walt
Subject: Re: how can I do Principal Components Analysis with octave?
Date: Wed, 12 May 2004 10:23:58 +0200
User-agent: Mozilla Thunderbird 0.5 (X11/20040306)

Maybe the attached klt.m will do the trick!

Regards
Stéfan

Henry F. Mollet wrote:
I believe you would need the programs doing the work. They are called
functions in Octave or Matlab, something like "PCA.m". I would be interested
myself so I checked using Google what MATLAB has:

Here is a list of the functions with a short description of each:
PRINCOMP - principal components from raw data matrix
PCACOV - pca from covariance matrix
PCARES - residuals from pca
BARTTEST - Bartlett's test for dimensionality.

Next I checked for the first two, namely "PRINCOMP" and "PCACOV" in
octave-forge but apparently neither is present. I guess we're out of luck
for the time being unless we have the capability to write the program.
Henry


on 5/11/04 1:04 AM, rino mailing at address@hidden wrote:


I'd like to do Principal Components Analysis with octave

What are the command I ave to write?

How to plot the result?



Thank you in advance for the time you spend to answer me, Mario.
##  KLT Perform the Karhunen-Loeve (a.k.a. Hotelling or PCA) transform.
##
##    [L, V] = klt(M, d or V[, options]);
##
##    M is the feature matrix (one feature per row) to be transformed.
##    The resulting matrix L has d dimensions (columns) unless
##    V (the transformation matrix) is given, in which case it has
##    the same nr of columns as V.
##
##    The transformation is defined as
##
##      y = V' * x;
##
##    or
##
##      y = V' * (x - mean_x)    (if option 'submean' is specified)
##
##    where x and y are column vectors and V the transformation
##    matrix.  V consists of the d largest eigenvectors.
##
##    The transform rotates data so that it lies along the
##    axis of greatest variance.  It is often used to
##    reduce dimensionality.
##
##    By default the eigenvalues are only plotted if no output
##    arguments are specified.  Use option 'plot' to override.
##
##    See Gonzales & Woods, Digital Image Processing
##    (p.675, Use of Principle Components)

## Author: Stefan van der Walt <address@hidden>, 2004

function [L, V] = klt(varargin)

    if (nargin < 2)
        usage("[L, V] = klt(M, d or V[, submean])");
    endif

    M = varargin{1};
    d = varargin{2};

    if (!ismatrix(M))
        error("klt: M must be a matrix");
    endif

    sub_mean = 0;
    plot_always = 0;
    if (length(varargin) > 2)
        for i = 1:length(varargin)
            if strcmpi(varargin{i}, "submean")
                sub_mean = 1;
            endif
            if strcmpi(varargin{i}, "plot")
                plot_always = 1;
            endif
        endfor
    endif

    if (!isscalar(d)) # V specified
        V = d;
        if (rows(d) != columns(M))
            error("klt: invalid dimensions for transformation matrix V (should 
be %d x D)", columns(M));
        endif
    else
        if (d > columns(M))
            warning("klt: asked for %d dimensions, but input only has %d!", [d 
columns(M)]);
            d = columns(M);
        endif
        
        C = cov(M);
    
        [V, lambda] = eig(C);
        [lambda, idx] = sort( diag(lambda) ); # sorted eigenvalues
        V = V(:, idx ); # sorted eigenvectors

        ## decreasing order
        lambda = flipud(lambda(:));
        V = fliplr(V);

        if ( (plot_always) || (nargout == 0) )      
            title("Eigenvalue decay");
            plot(lambda, ";Lambda;");
        endif
    
        ## choose the d eigenvectors corresponding to the d
        ## largest eigenvalues
        V = V(:, 1:d);
    endif

    ## transform every feature vector
    if ( !sub_mean )
        L = M * V;
    else
        L = ( M - repmat(mean(M), rows(M), 1) ) * V;    
    endif
    
endfunction

%!shared w

%!test
%!  [x, w] = klt([sqrt(2) sqrt(2); -2*sqrt(2) -2*sqrt(2)], 1);
%!  assert(x = [2; -4]);

%!test
%!  x = klt([sqrt(2) sqrt(2); -2*sqrt(2) -2*sqrt(2)], w);
%!  assert(x = [2; -2]);

%!test
%!  x = klt([sqrt(2) sqrt(2); -2*sqrt(2) -2*sqrt(2)], w, 'submean');
%!  assert(x = [3; -3]);

%!fail("klt(1)", "usage")

reply via email to

[Prev in Thread] Current Thread [Next in Thread]