octave_packages/nan-2.5.5/nantest.m

   1 % NANTEST checks several mathematical operations and a few
   2 % statistical functions for their correctness related to NaN's.
   3 % e.g. it checks norminv, normcdf, normpdf, sort, matrix division and multiplication.
   4 %
   5 %
   6 % see also: NANINSTTEST
   7 %
   8 % REFERENCE(S):
   9 % [1] W. Kahan (1996) Lecture notes on the Status of "IEEE Standard 754 for
  10 %     Binary Floating-point Arithmetic.
  11 %
  12
  13
  14 %    This program is free software; you can redistribute it and/or modify
  15 %    it under the terms of the GNU General Public License as published by
  16 %    the Free Software Foundation; either version 2 of the License, or
  17 %    (at your option) any later version.
  18 %
  19 %    This program is distributed in the hope that it will be useful,
  20 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 %    GNU General Public License for more details.
  23 %
  24 %    You should have received a copy of the GNU General Public License
  25 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  26
  27 %       $Id: nantest.m 8223 2011-04-20 09:16:06Z schloegl $
  28 %       Copyright (C) 2000-2004,2009 by Alois Schloegl <alois.schloegl@gmail.com>
  29 %       This script is part of the NaN-toolbox
  30 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  31
  32 %FLAG_WARNING = warning;
  33 %warning('off');
  34
  35 try
  36         x = randn([3,4,5]);
  37         x(~isnan(x)) = 0;
  38 catch
  39         fprintf(1,'WARNING: NANTEST fails for 3-DIM matrices. \n');
  40 end;
  41 try
  42         [s,n] = sumskipnan([nan,1,4,5]);
  43 catch
  44         fprintf(1,'WARNING: SUMSKIPNAN is not avaible. \n');
  45 end;
  46
  47 % check NORMPDF, NORMCDF, NORMINV
  48 x = [-inf,-2,-1,-.5,0,.5,1,2,3,inf,nan]';
  49 if exist('normpdf','file')==2,
  50         q(1) = sum(isnan(normpdf(x,2,0)))>sum(isnan(x));
  51         if q(1),
  52                 fprintf(1,'NORMPDF cannot handle v=0.\n');
  53                 fprintf(1,'-> NORMPDF should be replaced\n');
  54         end;
  55 end;
  56
  57 if exist('normcdf','file')==2,
  58         q(2) = sum(isnan(normcdf(x,2,0)))>sum(isnan(x));
  59         if q(2),
  60                 fprintf(1,'NORMCDF cannot handle v=0.\n');
  61                 fprintf(1,'-> NORMCDF should be replaced\n');
  62         end;
  63 end;
  64
  65 if exist('norminv','file')==2,
  66         p = [-inf,-.2,0,.2,.5,1,2,inf,nan];
  67         q(3) = sum(~isnan(norminv(p,2,0)))<4;
  68         if q(3),
  69                 fprintf(1,'NORMINV cannot handle  correctly v=0.\n');
  70                 fprintf(1,'-> NORMINV should be replaced\n');
  71         end;
  72         q(4) = ~isnan(norminv(0,NaN,0));
  73         q(5) = any(norminv(0.5,[1 2 3],0)~=(1:3));
  74 end;
  75
  76 if exist('tpdf','file')==2,
  77         q(6) = ~isnan(tpdf(nan,4));
  78         if q(6),
  79                 fprintf(1,'TPDF(NaN,4) does not return NaN\n');
  80                 fprintf(1,'-> TPDF should be replaced\n');
  81         end;
  82 end;
  83
  84 if exist('tcdf','file')==2,
  85         try
  86                 q(7) = ~isnan(tcdf(nan,4));
  87         catch
  88                 q(7) = 1;
  89         end;
  90         if q(7),
  91                 fprintf(1,'TCDF(NaN,4) does not return NaN\n');
  92                 fprintf(1,'-> TCDF should be replaced\n');
  93         end;
  94 end;
  95
  96 if exist('tinv','file')==2,
  97         try
  98                 q(8) = ~isnan(tinv(nan,4));
  99         catch
 100                 q(8) = 1;
 101         end;
 102         if q(8),
 103                 fprintf(1,'TINV(NaN,4) does not return NaN\n');
 104                 fprintf(1,'-> TINV should be replaced\n');
 105         end;
 106 end;
 107
 108 q(9) = isreal(double(2+3i));
 109 if q(9)
 110         printf('DOUBLE rejects imaginary part\n-> this can affect SUMSKIPNAN\n');
 111 end;
 112
 113 try
 114         x = reshape(1:6,3,2);
 115         [cc,nn] = covm(x+i*x,'e');
 116         q(10) = 0;
 117 catch
 118         q(10) = 1;
 119 end;
 120
 121 if 0,
 122 %%%%% MOD
 123 if exist('mod')>1,
 124         if (mod(5,0))~=0,
 125                 fprintf(1,'WARNING: MOD(x,0) does not return 0.\n');
 126         end;
 127         if isnan(mod(5,0)),
 128                 fprintf(1,'WARNING: MOD(x,0) returns NaN.\n');
 129         end;
 130         if isnan(mod(5,inf)),
 131                 fprintf(1,'WARNING: MOD(x,INF) returns NaN.\n');
 132         end;
 133 end;
 134 %%%%% REM
 135 if exist('rem')>1,
 136         if (rem(5,0))~=0,
 137                 fprintf(1,'WARNING: REM(x,0) does not return 0.\n');
 138         end;
 139         if isnan(rem(5,0)),
 140                 fprintf(1,'WARNING: REM(x,0) returns NaN.\n');
 141         end;
 142         if isnan(mod(5,inf)),
 143                 fprintf(1,'WARNING: REM(x,INF) returns NaN.\n');
 144         end;
 145 end;
 146 end;
 147
 148 %%%%% NANSUM(NAN) - this test addresses a problem in Matlab 5.3, 6.1 & 6.5
 149 if exist('nansum','file'),
 150         if isnan(nansum(nan)),
 151                 fprintf(1,'Warning: NANSUM(NaN) returns NaN instead of 0\n');
 152                 fprintf(1,'-> NANSUM should be replaced\n');
 153         end;
 154 end;
 155 %%%%% NANSUM(NAN) - this test addresses a problem in Matlab 5.3, 6.1 & 6.5
 156 if exist('nanstd','file'),
 157         if ~isnan(nanstd(0)),
 158                 fprintf(1,'Warning: NANSTD(x) with isscalar(x) returns 0 instead of NaN\n');
 159                 fprintf(1,'-> NANSTD should be replaced\n');
 160         end;
 161 end;
 162 %%%%% GEOMEAN - this test addresses a problem in Octave
 163 if exist('geomean','file'),
 164         if isnan(geomean((0:3)')),
 165                 fprintf(1,'Warning: GEOMEAN([0,1,2,3]) NaN instead of 0\n');
 166                 fprintf(1,'-> GEOMEAN should be replaced\n');
 167         end;
 168 end;
 169 %%%%% HARMMEAN - this test addresses a problem in Octave
 170 if exist('harmmean','file'),
 171         if isnan(harmmean(0:3)),
 172                 fprintf(1,'Warning: HARMMEAN([0,1,2,3]) NaN instead of 0\n');
 173                 fprintf(1,'-> HARMMEAN should be replaced\n');
 174         end;
 175 end;
 176 %%%%% BITAND - this test addresses a problem in Octave
 177 if exist('bitand')>1,
 178         if isnan(bitand(2^33-1,13)),
 179                 fprintf(1,'BITAND can return NaN. \n');
 180         end;
 181 end;
 182 %%%%% BITSHIFT - this test addresses a problem in Octave
 183 if exist('bitshift','file'),
 184         if isnan(bitshift(5,30,32)),
 185                 fprintf(1,'BITSHIFT can return NaN.\n');
 186         end;
 187 end;
 188 %%%%% ALL - this test addresses a problem in some old Octave and FreeMat v3.5
 189 if any(NaN)==1,
 190         fprintf(1,'WARNING: ANY(NaN) returns 1 instead of 0\n');
 191 end;
 192 if any([])==1,
 193         fprintf(1,'WARNING: ANY([]) returns 1 instead of 0\n');
 194 end;
 195 %%%%% ALL - this test addresses a problem in some old Octave and FreeMat v3.5
 196 if all(NaN)==0,
 197         fprintf(1,'WARNING: ALL(NaN) returns 0 instead of 1\n');
 198 end;
 199 if all([])==0,
 200         fprintf(1,'WARNING: ALL([]) returns 0 instead of 1\n');
 201 end;
 202
 203 %%%%% SORT - this was once a problem in Octave Version < 2.1.36 %%%%
 204 if ~all(isnan(sort([3,4,NaN,3,4,NaN]))==[0,0,0,0,1,1]),
 205         warning('Warning: SORT does not handle NaN.');
 206 end;
 207
 208 %%%%% commutativity of 0*NaN    %%% This test adresses a problem in Octave
 209 x=[-2:2;4:8]';
 210 y=x;y(2,1)=nan;y(4,2)=nan;
 211 B=[1,0,2;0,3,1];
 212 if ~all(all(isnan(y*B)==isnan(B'*y')')),
 213         fprintf(2,'WARNING: 0*NaN within matrix multiplication is not commutative\n');
 214 end;
 215
 216 % from Kahan (1996)
 217 tmp = (0-3*i)/inf;
 218 if isnan(tmp)
 219         fprintf(2,'WARNING: (0-3*i)/inf results in NaN instead of 0.\n');
 220 end;
 221
 222 %(roots([5,0,0])-[0;0])
 223 %(roots([2,-10,12])-[3;2])
 224 %(roots([2e-37,-2,2])-[1e37;1])
 225 %%%%% check nan/nan   %% this test addresses a problem in Matlab 5.3, 6.1 & 6.5
 226 p    = 4;
 227 tmp1 = repmat(nan,p)/repmat(nan,p);
 228 tmp2 = repmat(nan,p)\repmat(nan,p);
 229 tmp3 = repmat(0,p)/repmat(0,p);
 230 tmp4 = repmat(0,p)\repmat(0,p);
 231 tmp5 = repmat(0,p)*repmat(inf,p);
 232 tmp6 = repmat(inf,p)*repmat(0,p);
 233 x    = randn(100,1)*ones(1,p); y=x'*x;
 234 tmp7 = y/y;
 235 tmp8 = y\y;
 236
 237 if ~all(isnan(tmp1(:))),
 238         fprintf(1,'WARNING: matrix division NaN/NaN does not result in NaN\n');
 239 end;
 240 if ~all(isnan(tmp2(:))),
 241         fprintf(1,'WARNING: matrix division NaN\\NaN does not result in NaN\n');
 242 end;
 243 if ~all(isnan(tmp3(:))),
 244         fprintf(2,'WARNING: matrix division 0/0 does not result in NaN\n');
 245 end;
 246 if ~all(isnan(tmp4(:))),
 247         fprintf(2,'WARNING: matrix division 0\\0 does not result in NaN\n');
 248 end;
 249 if ~all(isnan(tmp5(:))),
 250         fprintf(2,'WARNING: matrix multiplication 0*inf does not result in NaN\n');
 251 end;
 252 if ~all(isnan(tmp6(:))),
 253         fprintf(2,'WARNING: matrix multiplication inf*0 does not result in NaN\n');
 254 end;
 255 if any(any(tmp7==inf));
 256         fprintf(2,'WARNING: right division of two singulare matrices return INF\n');
 257 end;
 258 if any(any(tmp8==inf));
 259         fprintf(2,'WARNING: left division of two singulare matrices return INF\n');
 260 end;
 261
 262 tmp  = [tmp1;tmp2;tmp3;tmp4;tmp5;tmp6;tmp7;tmp8];
 263
 264
 265
 266 %warning(FLAG_WARNING);
 267
 268
 269 %%%%% QUANTILE TEST
 270 d = [1 1 2 2 4 4 10 700];
 271 q = [-1,0,.05,.1,.25,.49,.5,.51,.75,.8, .999999,1,2];
 272 r = [ NaN, 1, 1, 1, 1.5, 2, 3, 4, 7, 10, 700, 700, NaN];
 273 if any( quantile(d, q) -  r>0)
 274         fprintf(1,'Quantile(1): failed\n');
 275 else
 276         fprintf(1,'Quantile(1): OK\n');
 277 end;
 278 if exist('histo3','file')
 279         H = histo3(d');
 280 else
 281         H.X = [1;2;4;10;700];
 282         H.H = [2;2;2;1;1];
 283         H.datatype = 'HISTOGRAM';
 284 end;
 285 if any( quantile(H, q)' -  r>0)
 286         fprintf(1,'Quantile(2): failed\n');
 287 else
 288         fprintf(1,'Quantile(2): OK\n');
 289 end;
 290