# cosmo corrΒΆ

```function c=cosmo_corr(x,y,corr_type)
% Computes correlation - faster than than matlab's "corr" for Pearson.
%
% c=comso_corr(x[,y[,corr_type]])
%
% Inputs:
%   x          PxM matrix.
%   y          PxN matrix (optional). If omitted then y=x.
%   corr_type  'Pearson' or 'Spearman' or 'Kendall' (optional). If omitted
%              then corrtype='Pearson' and the computation time is
%              significantly reduced for small matrices x and y (with
%              /tiny/ numerical imprecisions) by the use of a custom
%              implementation.
%              Using 'Kendall' required the matlab stats
%              toolbox and is currently not supported in Octave.
% Output:
%   c          MxN matrix with c(i,j)=corr(x(:,i),y(:,j),'type',corr_type).
%
% Notes:
%  - this function does not compute probability values.
%  - Using 'Kendall' for corr_type requires the matlab stats toolbox.
%
% Example:
%   % generate some pseudo-random data.
%   x=reshape(mod(2:7:100,41),5,[]);
%   y=reshape(mod(1:7:100,37),5,[]);
%   % compute builtin corr with cosmo_corr
%   % call the function first to avoid lookup delays; then measure time
%   c=corr(x,y);
%   cc=cosmo_corr(x,y);
%   % compute differences in output
%   delta=c-cc;
%   max_delta=max(abs(delta(:)));
%   fprintf('difference not greater than eps: %d\n',max_delta<=eps);
%   > difference not greater than eps: 1
%
%
% #   see the COPYING file distributed with CoSMoMVPA.           #

y_as_x=false;

if nargin<2
corr_type='Pearson';
y_as_x=true;
elseif ischar(y)
corr_type=y;
y_as_x=true;
elseif nargin<3
corr_type='Pearson';
end

if y_as_x
y=x;
end

switch corr_type
case 'Pearson'
c=corr_pearson(x,y,y_as_x);

case 'Spearman'
x_ranks=cosmo_tiedrank(x);
y_ranks=cosmo_tiedrank(y);

c=corr_pearson(x_ranks,y_ranks,y_as_x);

otherwise
% fall-back: use Matlab's function
% will puke if no Matlab stat toolbox
c=corr(x,y,'type',corr_type);
end

function c=corr_pearson(x,y,y_as_x)

% speed-optimized version
nx=size(x,1);
ny=size(y,1);

% subtract mean
xd=bsxfun(@minus,x,sum(x,1)/nx);
yd=bsxfun(@minus,y,sum(y,1)/ny);

% normalization
n=1/(size(x,1)-1);

% standard deviation
xs=(n*sum(xd .^ 2)).^-0.5;
ys=(n*sum(yd .^ 2)).^-0.5;

% compute correlations
c=n * (xd' * yd) .* (xs' * ys);

if y_as_x
% ensure diagonal elements are 1
c=(c+c')*.5;
dc=diag(c);
c=(c-diag(dc))+eye(numel(dc));
end
```