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