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