octave_packages/nan-2.5.5/cor.m

   1 function [r2] = cor(X,Y);
   2 % COR calculates the correlation matrix
   3 % X and Y can contain missing values encoded with NaN.
   4 % NaN's are skipped, NaN do not result in a NaN output.
   5 % (Its assumed that the occurence of NaN's is uncorrelated)
   6 % The output gives NaN only if there are insufficient input data
   7 %
   8 % COR(X);
   9 %      calculates the (auto-)correlation matrix of X
  10 % COR(X,Y);
  11 %      calculates the crosscorrelation between X and Y
  12 %
  13 % c = COR(...);
  14 %       c is the correlation matrix
  15 %
  16 % W     weights to compute weighted mean (default: [])
  17 %       if W=[], all weights are 1.
  18 %       number of elements in W must match size(x,DIM)
  19
  20 % NOTE: Under certain circumstances (Missing values and small number of samples)
  21 %   abs(COR) can be larger than 1.
  22 %   If you need abs(COR)<=1, use CORRCOEF. CORRCOEF garantees abs(COR)<=1.
  23 %
  24 % see also: SUMSKIPNAN, COVM, COV, CORRCOEF
  25 %
  26 % REFERENCES:
  27 % http://mathworld.wolfram.com/CorrelationCoefficient.html
  28
  29
  30 %       $Id: cor.m 8223 2011-04-20 09:16:06Z schloegl $
  31 %       Copyright (C) 2000-2004,2010 by Alois Schloegl <alois.schloegl@gmail.com>
  32 %       This function is part of the NaN-toolbox
  33 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  34
  35
  36 %    This program is free software; you can redistribute it and/or modify
  37 %    it under the terms of the GNU General Public License as published by
  38 %    the Free Software Foundation; either version 2 of the License, or
  39 %    (at your option) any later version.
  40 %
  41 %    This program is distributed in the hope that it will be useful,
  42 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  43 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  44 %    GNU General Public License for more details.
  45 %
  46 %    You should have received a copy of the GNU General Public License
  47 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  48
  49
  50 if nargin==1
  51         Y = [];
  52 elseif nargin==0
  53         fprintf(2,'Error COR: Missing argument(s)\n');
  54 end;
  55
  56 [r1,c1]=size(X);
  57 if (c1>r1),
  58         fprintf(2,'Warning COR: Covariance is ill-defined, because of too less observations (rows).\n');
  59 end;
  60
  61 [r1,c1]=size(X);
  62 if ~isempty(Y)
  63         [r2,c2]=size(Y);
  64         if r1~=r2,
  65                 fprintf(2,'Error COR: X and Y must have the same number of observations (rows).\n');
  66                 return;
  67         end;
  68 else
  69         [r2,c2]=size(X);
  70 end;
  71
  72 if (c1>r1) || (c2>r2),
  73         fprintf(2,'Warning COR: Covariance is ill-defined, because of too less observations (rows).\n');
  74 end;
  75
  76 if ~isempty(Y),
  77         [S1,N1,SSQ1] = sumskipnan(X,1);
  78         [S2,N2,SSQ2] = sumskipnan(Y,1);
  79
  80         NN = double(~isnan(X)')*double(~isnan(Y));
  81         X(isnan(X)) = 0; % skip NaN's
  82         Y(isnan(Y)) = 0; % skip NaN's
  83         CC = X'*Y;
  84
  85         M1 = S1./N1;
  86         M2 = S2./N2;
  87         cc = CC./NN - M1'*M2;
  88         r2 = cc./sqrt((SSQ1./N1-M1.*M1)'*(SSQ2./N2-M2.*M2));
  89
  90 else
  91         [S,N,SSQ] = sumskipnan(X,1);
  92
  93         NN = double(~isnan(X)')*double(~isnan(X));
  94         X(isnan(X)) = 0; % skip NaN's
  95         CC = X'*X;
  96
  97         M  = S./N;
  98         cc = CC./NN - M'*M;
  99         v  = (SSQ./N- M.*M);  %max(N-1,0);
 100         r2 = cc./sqrt(v'*v);
 101 end;