Mercurial > octave-nkf

comparison scripts/statistics/base/skewness.m @ 17586:8334144458a4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
skewness.m: Fix "division by zero" warnings. Improve input validation for FLAG. * scripts/statistics/base/skewness.m: Only divide by std^3 when it is non-zero. Check that FLAG is scalar before using. Increase test tolerance to 2*eps to avoid failure on some platforms. Add %!test for division-by-zero warnings.
author Julien Bect <julien.bect@supelec.fr>
date Mon, 07 Oct 2013 15:21:06 +0200
parents 4cb05034f1c6
children a13ff4521538
comparison
equal deleted inserted replaced
17585:36b9fa789d8e 17586:8334144458a4
89 endif 89 endif
90 90
91 if (nargin < 2 || isempty (flag)) 91 if (nargin < 2 || isempty (flag))
92 flag = 1; # default: do not use the "bias corrected" version 92 flag = 1; # default: do not use the "bias corrected" version
93 else 93 else
94 flag = double (flag); 94 if ((! isscalar (flag)) || (flag != 0 && flag != 1))
95 if (flag != 0 && flag != 1)
96 error ("skewness: FLAG must be 0 or 1"); 95 error ("skewness: FLAG must be 0 or 1");
97 endif 96 endif
98 endif 97 endif
99 98
100 nd = ndims (x); 99 nd = ndims (x);
112 sz(dim) = 1; 111 sz(dim) = 1;
113 112
114 x = center (x, dim); # center also promotes integer, logical to double 113 x = center (x, dim); # center also promotes integer, logical to double
115 s = std (x, 1, dim); # Normalize with 1/N 114 s = std (x, 1, dim); # Normalize with 1/N
116 y = sum (x .^ 3, dim); 115 y = sum (x .^ 3, dim);
117 y ./= (n * s .^ 3); 116 idx = (s != 0);
118 y(s == 0) = NaN; 117 y(idx) ./= (n * s(idx) .^ 3);
118 y(! idx) = NaN;
119 119
120 ## Apply bias correction to the third central sample moment 120 ## Apply bias correction to the third central sample moment
121 if (flag == 0) 121 if (flag == 0)
122 if (n > 2) 122 if (n > 2)
123 y *= sqrt (n * (n - 1)) / (n - 2); 123 y *= sqrt (n * (n - 1)) / (n - 2);
139 %! x = [0; 0; 0; 1]; 139 %! x = [0; 0; 0; 1];
140 %! y = [x, 2*x]; 140 %! y = [x, 2*x];
141 %! assert (skewness (y), 1.154700538379251 * [1 1], 5*eps); 141 %! assert (skewness (y), 1.154700538379251 * [1 1], 5*eps);
142 142
143 %!assert (skewness ([1:5 10; 1:5 10], 0, 2), 1.439590274527954 * [1; 1], eps) 143 %!assert (skewness ([1:5 10; 1:5 10], 0, 2), 1.439590274527954 * [1; 1], eps)
144 %!assert (skewness ([1:5 10; 1:5 10], 1, 2), 1.051328089232020 * [1; 1], eps) 144 %!assert (skewness ([1:5 10; 1:5 10], 1, 2), 1.051328089232020 * [1; 1], 2*eps)
145 %!assert (skewness ([1:5 10; 1:5 10], [], 2), 1.051328089232020 * [1; 1], eps) 145 %!assert (skewness ([1:5 10; 1:5 10], [], 2), 1.051328089232020 * [1; 1], 2*eps)
146 146
147 ## Test behaviour on single input 147 ## Test behaviour on single input
148 %!assert (skewness (single ([1:5 10])), single (1.0513283), eps ("single")) 148 %!assert (skewness (single ([1:5 10])), single (1.0513283), eps ("single"))
149 %!assert (skewness (single ([1 2]), 0), single (NaN)) 149 %!assert (skewness (single ([1 2]), 0), single (NaN))
150
151 ## Verify no "division-by-zero" warnings
152 %!test
153 %! wstate = warning ("query", "Octave:divide-by-zero");
154 %! unwind_protect
155 %! lastwarn (""); # clear last warning
156 %! skewness (1);
157 %! assert (lastwarn (), "");
158 %! unwind_protect_cleanup
159 %! warning (wstate, "Octave:divide-by-zero");
160 %! end_unwind_protect
150 161
151 ## Test input validation 162 ## Test input validation
152 %!error skewness () 163 %!error skewness ()
153 %!error skewness (1, 2, 3) 164 %!error skewness (1, 2, 3)
154 %!error <X must be a numeric vector or matrix> skewness (['A'; 'B']) 165 %!error <X must be a numeric vector or matrix> skewness (['A'; 'B'])
155 %!error <FLAG must be 0 or 1> skewness (1, 2) 166 %!error <FLAG must be 0 or 1> skewness (1, 2)
167 %!error <FLAG must be 0 or 1> skewness (1, [1 0])
156 %!error <DIM must be an integer> skewness (1, [], ones (2,2)) 168 %!error <DIM must be an integer> skewness (1, [], ones (2,2))
157 %!error <DIM must be an integer> skewness (1, [], 1.5) 169 %!error <DIM must be an integer> skewness (1, [], 1.5)
158 %!error <DIM must be .* a valid dimension> skewness (1, [], 0) 170 %!error <DIM must be .* a valid dimension> skewness (1, [], 0)
159 %!error <DIM must be .* a valid dimension> skewness (1, [], 3) 171 %!error <DIM must be .* a valid dimension> skewness (1, [], 3)
160 172