Mercurial > octave-nkf

comparison scripts/statistics/base/skewness.m @ 17582:7004c733412f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
skewness.m: Update cset 3ef7d28833f3 to Octave coding conventions. * scripts/statistics/base/skewness.m: Use @group for multi-line examples in docstring. Use Matlab-compatible definition of bias correction factor which is "sqrt (N*(N-1))/(N-2)". Also accept logical vectors. Use parentheses around conditions in if statements. Use '!' instead of '~' for logical negation. Use in-place operators for efficiency. Use 'NaN' instead of 'nan'. Use eps instead of hardcoded tolerance. Look for error messages in %!error blocks.
author Rik <rik@octave.org>
date Sun, 06 Oct 2013 16:16:36 -0700
parents ab75e72c5b36
children 4cb05034f1c6
comparison
equal deleted inserted replaced
17581:ab75e72c5b36 17582:7004c733412f
31 ## its (uncorrected) standard deviation. 31 ## its (uncorrected) standard deviation.
32 ## @end tex 32 ## @end tex
33 ## @ifnottex 33 ## @ifnottex
34 ## 34 ##
35 ## @example 35 ## @example
36 ## @group
36 ## mean ((@var{x} - mean (@var{x})).^3) 37 ## mean ((@var{x} - mean (@var{x})).^3)
37 ## skewness (@var{X}) = ------------------------. 38 ## skewness (@var{X}) = ------------------------.
38 ## std (@var{x}, 1).^3 39 ## std (@var{x}).^3
40 ## @end group
39 ## @end example 41 ## @end example
40 ## 42 ##
41 ## @end ifnottex 43 ## @end ifnottex
42 ## 44 ##
43 ## @noindent 45 ## @noindent
44 ## The optional argument @var{flag} controls which normalization is used. 46 ## The optional argument @var{flag} controls which normalization is used.
45 ## If @var{flag} is equal to 1 (default value, used when @var{flag} is omitted 47 ## If @var{flag} is equal to 1 (default value, used when @var{flag} is omitted
46 ## or empty), return the sample skewness as defined above. If @var{flag} is 48 ## or empty), return the sample skewness as defined above. If @var{flag} is
47 ## equal to 0, return the adjusted skewness coefficient instead: 49 ## equal to 0, return the adjusted skewness coefficient instead:
48 ## @tex 50 ## @tex
49 ## $$ 51 ## $$
50 ## {\rm skewness} (@var{x}) = {{{N \over (N - 1) (N - 2)}\, 52 ## {\rm skewness} (@var{x}) = {\sqrt{N (N - 1)} \over N - 2} \times \,
51 ## \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^3} \over @var{s}^3}, 53 ## {{{1 \over N} \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^3} \over \sigma^3}
52 ## $$ 54 ## $$
53 ## where @var{s} is the corrected standard deviation of @var{x}.
54 ## @end tex 55 ## @end tex
55 ## @ifnottex 56 ## @ifnottex
56 ## 57 ##
57 ## @example 58 ## @example
58 ## N^2 mean ((@var{x} - mean (@var{x})).^3) 59 ## @group
60 ## sqrt (N*(N-1)) mean ((@var{x} - mean (@var{x})).^3)
59 ## skewness (@var{X}, 0) = -------------- * ------------------------. 61 ## skewness (@var{X}, 0) = -------------- * ------------------------.
60 ## (N - 1)(N - 2) std (@var{x}, 0).^3 62 ## (N - 2) std (@var{x}).^3
63 ## @end group
61 ## @end example 64 ## @end example
62 ## 65 ##
63 ## @end ifnottex 66 ## @end ifnottex
64 ## The adjusted skewness coefficient is obtained by replacing the sample second 67 ## The adjusted skewness coefficient is obtained by replacing the sample second
65 ## and third central moments by their bias-corrected versions. 68 ## and third central moments by their bias-corrected versions.
66 ## 69 ##
67 ## If @var{x} is a matrix, or more generally a multidimensional array, return 70 ## If @var{x} is a matrix, or more generally a multi-dimensional array, return
68 ## the skewness along the first non-singleton dimension. If the optional 71 ## the skewness along the first non-singleton dimension. If the optional
69 ## @var{dim} argument is given, operate along this dimension. 72 ## @var{dim} argument is given, operate along this dimension.
70 ## 73 ##
71 ## @seealso{var, kurtosis, moment} 74 ## @seealso{var, kurtosis, moment}
72 ## @end deftypefn 75 ## @end deftypefn
79 82
80 if (nargin < 1) || (nargin > 3) 83 if (nargin < 1) || (nargin > 3)
81 print_usage (); 84 print_usage ();
82 endif 85 endif
83 86
84 if (! isnumeric (x)) || islogical (x) 87 if (! (isnumeric (x) || islogical (x)))
85 error ("skewness: X must be a numeric vector or matrix"); 88 error ("skewness: X must be a numeric vector or matrix");
86 endif 89 endif
87 90
88 if (nargin < 2) || isempty (flag) 91 if (nargin < 2 || isempty (flag))
89 flag = 1; # default: do not use the "bias corrected" version 92 flag = 1; # default: do not use the "bias corrected" version
90 else 93 else
91 flag = double (flag); 94 flag = double (flag);
92 if (~ isequal (flag, 0)) && (~ isequal (flag, 1)) 95 if (flag != 0 && flag != 1)
93 error ("skewness: FLAG must be 0 or 1"); 96 error ("skewness: FLAG must be 0 or 1");
94 end 97 end
95 endif 98 endif
96 99
97 nd = ndims (x); 100 nd = ndims (x);
98 sz = size (x); 101 sz = size (x);
99 if nargin < 3 102 if (nargin < 3)
100 ## Find the first non-singleton dimension. 103 ## Find the first non-singleton dimension.
101 (dim = find (sz > 1, 1)) || (dim = 1); 104 (dim = find (sz > 1, 1)) || (dim = 1);
102 else 105 else
103 if (!(isscalar (dim) && dim == fix (dim)) || !(1 <= dim && dim <= nd)) 106 if (!(isscalar (dim) && dim == fix (dim)) || !(1 <= dim && dim <= nd))
104 error ("skewness: DIM must be an integer and a valid dimension"); 107 error ("skewness: DIM must be an integer and a valid dimension");
106 endif 109 endif
107 110
108 n = sz(dim); 111 n = sz(dim);
109 sz(dim) = 1; 112 sz(dim) = 1;
110 113
111 ## In the following sequence of computations, if x is of class single, 114 x = center (x, dim); # center also promotes integer, logical to double
112 ## then the result y is also of class single 115 s = std (x, 1, dim); # Normalize with 1/N
113 x = center (x, dim);
114 s = std (x, flag, dim); # use the biased estimator of the variance
115 y = sum (x .^ 3, dim); 116 y = sum (x .^ 3, dim);
116 y = y ./ (n * s .^ 3); 117 y ./= (n * s .^ 3);
117 y(s <= 0) = nan; 118 y(s == 0) = NaN;
118 119
119 ## Apply bias correction to the third central sample moment 120 ## Apply bias correction to the third central sample moment
120 if flag == 0 121 if (flag == 0)
121 if n > 2 122 if (n > 2)
122 y = y * (n * n ) / ((n - 1) * (n - 2)); 123 y *= sqrt (n * (n - 1)) / (n - 2);
123 else 124 else
124 y(:) = nan; 125 y(:) = NaN;
125 end 126 end
126 endif 127 endif
127 128
128 endfunction 129 endfunction
129 130
130 131
131 %!assert (skewness ([-1, 0, 1]), 0) 132 %!assert (skewness ([-1, 0, 1]), 0)
132 %!assert (skewness ([-2, 0, 1]) < 0) 133 %!assert (skewness ([-2, 0, 1]) < 0)
133 %!assert (skewness ([-1, 0, 2]) > 0) 134 %!assert (skewness ([-1, 0, 2]) > 0)
134 %!assert (skewness ([-3, 0, 1]) == -1 * skewness ([-1, 0, 3])) 135 %!assert (skewness ([-3, 0, 1]) == -1 * skewness ([-1, 0, 3]))
135 %!assert (skewness (ones (3, 5)), nan (1, 5)) 136 %!assert (skewness (ones (3, 5)), NaN (1, 5))
136 137
137 %!test 138 %!test
138 %! x = [0; 0; 0; 1]; 139 %! x = [0; 0; 0; 1];
139 %! y = [x, 2*x]; 140 %! y = [x, 2*x];
140 %! assert (skewness (y), 1.154700538379251 * [1 1], 1e-13) 141 %! assert (skewness (y), 1.154700538379251 * [1 1], 5*eps);
141 142
142 %!assert (skewness ([1:5 10; 1:5 10], 0, 2), 1.439590274527954 * [1; 1], 1e-15) 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], 1, 2), 1.051328089232020 * [1; 1], 1e-15) 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], [], 2), 1.051328089232020 * [1; 1], 1e-15) 145 %!assert (skewness ([1:5 10; 1:5 10], [], 2), 1.051328089232020 * [1; 1], eps)
145 146
146 ## Test behaviour on single input 147 ## Test behaviour on single input
147 %! assert (skewness (single ([1:5 10])), single (1.0513283)) 148 %!assert (skewness (single ([1:5 10])), single (1.0513283), eps ("single"))
148 %! assert (skewness (single ([1 2]), 0), single (nan)) 149 %!assert (skewness (single ([1 2]), 0), single (NaN))
149 150
150 ## Test input validation 151 ## Test input validation
151 %!error skewness () 152 %!error skewness ()
152 %!error skewness (1, 2, 3) 153 %!error skewness (1, 2, 3)
153 %!error skewness (['A'; 'B']) 154 %!error <X must be a numeric vector or matrix> skewness (['A'; 'B'])
154 %!error skewness (1, ones (2,2)) 155 %!error <FLAG must be 0 or 1> skewness (1, 2)
155 %!error skewness (1, 1.5) 156 %!error <DIM must be an integer> skewness (1, [], ones (2,2))
156 %!error skewness (1, [], 0) 157 %!error <DIM must be an integer> skewness (1, [], 1.5)
157 %!error skewness (1, 3) 158 %!error <DIM must be .* a valid dimension> skewness (1, [], 0)
159 %!error <DIM must be .* a valid dimension> skewness (1, [], 3)
158 160