octave_packages/nan-2.5.5/covm.m

   1 function [CC,NN] = covm(X,Y,Mode,W)
   2 % COVM generates covariance 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 % The output gives NaN only if there are insufficient input data
   6 %
   7 % COVM(X,Mode);
   8 %      calculates the (auto-)correlation matrix of X
   9 % COVM(X,Y,Mode);
  10 %      calculates the crosscorrelation between X and Y
  11 % COVM(...,W);
  12 %       weighted crosscorrelation
  13 %
  14 % Mode = 'M' minimum or standard mode [default]
  15 %       C = X'*X; or X'*Y correlation matrix
  16 %
  17 % Mode = 'E' extended mode
  18 %       C = [1 X]'*[1 X]; % l is a matching column of 1's
  19 %       C is additive, i.e. it can be applied to subsequent blocks and summed up afterwards
  20 %       the mean (or sum) is stored on the 1st row and column of C
  21 %
  22 % Mode = 'D' or 'D0' detrended mode
  23 %       the mean of X (and Y) is removed. If combined with extended mode (Mode='DE'),
  24 %       the mean (or sum) is stored in the 1st row and column of C.
  25 %       The default scaling is factor (N-1).
  26 % Mode = 'D1' is the same as 'D' but uses N for scaling.
  27 %
  28 % C = covm(...);
  29 %       C is the scaled by N in Mode M and by (N-1) in mode D.
  30 % [C,N] = covm(...);
  31 %       C is not scaled, provides the scaling factor N
  32 %       C./N gives the scaled version.
  33 %
  34 % see also: DECOVM, XCOVF
  35
  36 %       $Id: covm.m 9032 2011-11-08 20:25:36Z schloegl $
  37 %       Copyright (C) 2000-2005,2009 by Alois Schloegl <alois.schloegl@gmail.com>
  38 %       This function is part of the NaN-toolbox
  39 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  40
  41 %    This program is free software; you can redistribute it and/or modify
  42 %    it under the terms of the GNU General Public License as published by
  43 %    the Free Software Foundation; either version 3 of the License, or
  44 %    (at your option) any later version.
  45 %
  46 %    This program is distributed in the hope that it will be useful,
  47 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  48 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  49 %    GNU General Public License for more details.
  50 %
  51 %    You should have received a copy of the GNU General Public License
  52 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  53
  54
  55 global FLAG_NANS_OCCURED;
  56
  57 if nargin<3,
  58         W = [];
  59         if nargin==2,
  60                 if isnumeric(Y),
  61                         Mode='M';
  62                 else
  63                         Mode=Y;
  64                         Y=[];
  65                 end;
  66         elseif nargin==1,
  67                 Mode = 'M';
  68                 Y = [];
  69         elseif nargin==0,
  70                 error('Missing argument(s)');
  71         end;
  72
  73 elseif (nargin==3) && isnumeric(Y) && ~isnumeric(Mode);
  74         W = [];
  75
  76 elseif (nargin==3) && ~isnumeric(Y) && isnumeric(Mode);
  77         W = Mode;
  78         Mode = Y;
  79         Y = [];
  80
  81 elseif (nargin==4) && ~isnumeric(Mode) && isnumeric(Y);
  82         ; %% ok
  83 else
  84         error('invalid input arguments');
  85 end;
  86
  87 Mode = upper(Mode);
  88
  89 [r1,c1]=size(X);
  90 if ~isempty(Y)
  91         [r2,c2]=size(Y);
  92         if r1~=r2,
  93                 error('X and Y must have the same number of observations (rows).');
  94         end;
  95 else
  96         [r2,c2]=size(X);
  97 end;
  98
  99 persistent mexFLAG2;
 100 persistent mexFLAG;
 101 if isempty(mexFLAG2)
 102         mexFLAG2 = exist('covm_mex','file');
 103 end;
 104 if isempty(mexFLAG)
 105         mexFLAG = exist('sumskipnan_mex','file');
 106 end;
 107
 108
 109 if ~isempty(W)
 110         W = W(:);
 111         if (r1~=numel(W))
 112                 error('Error COVM: size of weight vector does not fit number of rows');
 113         end;
 114         %w = spdiags(W(:),0,numel(W),numel(W));
 115         %nn = sum(W(:));
 116         nn = sum(W);
 117 else
 118         nn = r1;
 119 end;
 120
 121
 122 if mexFLAG2 && mexFLAG && ~isempty(W),
 123         %% the mex-functions here are much slower than the m-scripts below
 124         %% however, the mex-functions support weighting of samples.
 125         if isempty(FLAG_NANS_OCCURED),
 126                 %% mex-files require that FLAG_NANS_OCCURED is not empty,
 127                 %% otherwise, the status of NAN occurence can not be returned.
 128                 FLAG_NANS_OCCURED = logical(0);  % default value
 129         end;
 130
 131         if any(Mode=='D') || any(Mode=='E'),
 132                 [S1,N1] = sumskipnan(X,1,W);
 133                 if ~isempty(Y)
 134                         [S2,N2] = sumskipnan(Y,1,W);
 135                 else
 136                         S2 = S1; N2 = N1;
 137                 end;
 138                 if any(Mode=='D'), % detrending mode
 139                         X  = X - ones(r1,1)*(S1./N1);
 140                         if ~isempty(Y)
 141                                 Y  = Y - ones(r1,1)*(S2./N2);
 142                         end;
 143                 end;
 144         end;
 145
 146         [CC,NN] = covm_mex(real(X), real(Y), FLAG_NANS_OCCURED, W);
 147         %% complex matrices
 148         if ~isreal(X) && ~isreal(Y)
 149                 [iCC,inn] = covm_mex(imag(X), imag(Y), FLAG_NANS_OCCURED, W);
 150                 CC = CC + iCC;
 151         end;
 152         if isempty(Y) Y = X; end;
 153         if ~isreal(X)
 154                 [iCC,inn] = covm_mex(imag(X), real(Y), FLAG_NANS_OCCURED, W);
 155                 CC = CC - i*iCC;
 156         end;
 157         if ~isreal(Y)
 158                 [iCC,inn] = covm_mex(real(X), imag(Y), FLAG_NANS_OCCURED, W);
 159                 CC = CC + i*iCC;
 160         end;
 161
 162         if any(Mode=='D') && ~any(Mode=='1'),  %  'D1'
 163                 NN = max(NN-1,0);
 164         end;
 165         if any(Mode=='E'), % extended mode
 166                 NN = [nn, N2; N1', NN];
 167                 CC = [nn, S2; S1', CC];
 168         end;
 169
 170
 171 elseif ~isempty(W),
 172
 173         error('Error COVM: weighted COVM requires sumskipnan_mex and covm_mex but it is not available');
 174
 175         %% weighted covm without mex-file support
 176         %% this part is not working.
 177
 178 elseif ~isempty(Y),
 179         if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M
 180                 NN = real(X==X)'*real(Y==Y);
 181                 FLAG_NANS_OCCURED = any(NN(:)<nn);
 182                 X(X~=X) = 0; % skip NaN's
 183                 Y(Y~=Y) = 0; % skip NaN's
 184                 CC = X'*Y;
 185
 186         else  % if any(Mode=='D') | any(Mode=='E'),
 187                 [S1,N1] = sumskipnan(X,1);
 188                 [S2,N2] = sumskipnan(Y,1);
 189                 NN = real(X==X)'*real(Y==Y);
 190
 191                 if any(Mode=='D'), % detrending mode
 192                         X  = X - ones(r1,1)*(S1./N1);
 193                         Y  = Y - ones(r1,1)*(S2./N2);
 194                         if any(Mode=='1'),  %  'D1'
 195                                 NN = NN;
 196                         else    %  'D0'
 197                                 NN = max(NN-1,0);
 198                         end;
 199                 end;
 200                 X(X~=X) = 0; % skip NaN's
 201                 Y(Y~=Y) = 0; % skip NaN's
 202                 CC = X'*Y;
 203
 204                 if any(Mode=='E'), % extended mode
 205                         NN = [nn, N2; N1', NN];
 206                         CC = [nn, S2; S1', CC];
 207                 end;
 208         end;
 209
 210 else
 211         if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M
 212                 tmp = real(X==X);
 213                 NN  = tmp'*tmp;
 214                 X(X~=X) = 0; % skip NaN's
 215                 CC = X'*X;
 216                 FLAG_NANS_OCCURED = any(NN(:)<nn);
 217
 218         else  % if any(Mode=='D') | any(Mode=='E'),
 219                 [S,N] = sumskipnan(X,1);
 220                 tmp = real(X==X);
 221                 NN  = tmp'*tmp;
 222                 if any(Mode=='D'), % detrending mode
 223                         X  = X - ones(r1,1)*(S./N);
 224                         if any(Mode=='1'),  %  'D1'
 225                                 NN = NN;
 226                         else  %  'D0'
 227                                 NN = max(NN-1,0);
 228                         end;
 229                 end;
 230
 231                 X(X~=X) = 0; % skip NaN's
 232                 CC = X'*X;
 233                 if any(Mode=='E'), % extended mode
 234                         NN = [nn, N; N', NN];
 235                         CC = [nn, S; S', CC];
 236                 end;
 237         end
 238
 239 end;
 240
 241
 242 if nargout<2
 243         CC = CC./NN; % unbiased
 244 end;
 245
 246 %!assert(covm([1;NaN;2],'D'),0.5)
 247 %!assert(covm([1;NaN;2],'M'),2.5)
 248 %!assert(covm([1;NaN;2],'E'),[1,1.5;1.5,2.5])