Mercurial > octave-nkf

comparison scripts/statistics/base/skewness.m @ 17581:ab75e72c5b36

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
skewness.m: Improve compatibility with Matlab's skewness function * skewness.m: Add a 'flag' input argument to select either the sample skewness (default) or its "bias corrected" version. Return NaN (instead of 0) when the standard deviation is zero. Fix documentation. Fix old tests and add some new ones.
author Julien Bect <julien.bect@supelec.fr>
date Tue, 01 Oct 2013 21:51:17 +0200
parents f3d52523cde1
children 7004c733412f
comparison
equal deleted inserted replaced
17580:3ef7d28833f3 17581:ab75e72c5b36
1 ## Copyright (C) 2013 Julien Bect
1 ## Copyright (C) 1996-2012 John W. Eaton 2 ## Copyright (C) 1996-2012 John W. Eaton
2 ## 3 ##
3 ## This file is part of Octave. 4 ## This file is part of Octave.
4 ## 5 ##
5 ## Octave is free software; you can redistribute it and/or modify it 6 ## Octave is free software; you can redistribute it and/or modify it
16 ## along with Octave; see the file COPYING. If not, see 17 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 18 ## <http://www.gnu.org/licenses/>.
18 19
19 ## -*- texinfo -*- 20 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {} skewness (@var{x}) 21 ## @deftypefn {Function File} {} skewness (@var{x})
21 ## @deftypefnx {Function File} {} skewness (@var{x}, @var{dim}) 22 ## @deftypefnx {Function File} {} skewness (@var{x}, @var{flag})
22 ## Compute the skewness of the elements of the vector @var{x}. 23 ## @deftypefnx {Function File} {} skewness (@var{x}, @var{flag}, @var{dim})
24 ## Compute the sample skewness of the elements of @var{x}:
23 ## @tex 25 ## @tex
24 ## $$ 26 ## $$
25 ## {\rm skewness} (x) = {1\over N \sigma^3} \sum_{i=1}^N (x_i-\bar{x})^3 27 ## {\rm skewness} (@var{x}) = {{{1\over N}\,
28 ## \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^3} \over \sigma^3},
26 ## $$ 29 ## $$
27 ## where $\bar{x}$ is the mean value of $x$. 30 ## where $N$ is the length of @var{x}, $\bar{@var{x}}$ its mean and $\sigma$
31 ## its (uncorrected) standard deviation.
28 ## @end tex 32 ## @end tex
29 ## @ifnottex 33 ## @ifnottex
30 ## 34 ##
31 ## @example 35 ## @example
32 ## skewness (x) = 1/N std(x)^(-3) sum ((x - mean(x)).^3) 36 ## mean ((@var{x} - mean (@var{x})).^3)
37 ## skewness (@var{X}) = ------------------------.
38 ## std (@var{x}, 1).^3
33 ## @end example 39 ## @end example
34 ## 40 ##
35 ## @end ifnottex 41 ## @end ifnottex
36 ## 42 ##
37 ## @noindent 43 ## @noindent
38 ## If @var{x} is a matrix, return the skewness along the 44 ## The optional argument @var{flag} controls which normalization is used.
39 ## first non-singleton dimension of the matrix. If the optional 45 ## 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
47 ## equal to 0, return the adjusted skewness coefficient instead:
48 ## @tex
49 ## $$
50 ## {\rm skewness} (@var{x}) = {{{N \over (N - 1) (N - 2)}\,
51 ## \sum_{i=1}^N (@var{x}_i - \bar{@var{x}})^3} \over @var{s}^3},
52 ## $$
53 ## where @var{s} is the corrected standard deviation of @var{x}.
54 ## @end tex
55 ## @ifnottex
56 ##
57 ## @example
58 ## N^2 mean ((@var{x} - mean (@var{x})).^3)
59 ## skewness (@var{X}, 0) = -------------- * ------------------------.
60 ## (N - 1)(N - 2) std (@var{x}, 0).^3
61 ## @end example
62 ##
63 ## @end ifnottex
64 ## The adjusted skewness coefficient is obtained by replacing the sample second
65 ## and third central moments by their bias-corrected versions.
66 ##
67 ## If @var{x} is a matrix, or more generally a multidimensional array, return
68 ## the skewness along the first non-singleton dimension. If the optional
40 ## @var{dim} argument is given, operate along this dimension. 69 ## @var{dim} argument is given, operate along this dimension.
70 ##
41 ## @seealso{var, kurtosis, moment} 71 ## @seealso{var, kurtosis, moment}
42 ## @end deftypefn 72 ## @end deftypefn
43 73
44 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> 74 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
45 ## Created: 29 July 1994 75 ## Created: 29 July 1994
46 ## Adapted-By: jwe 76 ## Adapted-By: jwe
47 77
48 function retval = skewness (x, dim) 78 function y = skewness (x, flag, dim)
49 79
50 if (nargin != 1 && nargin != 2) 80 if (nargin < 1) || (nargin > 3)
51 print_usage (); 81 print_usage ();
52 endif 82 endif
53 83
54 if (! (isnumeric (x) || islogical (x))) 84 if (! isnumeric (x)) || islogical (x)
55 error ("skewness: X must be a numeric vector or matrix"); 85 error ("skewness: X must be a numeric vector or matrix");
86 endif
87
88 if (nargin < 2) || isempty (flag)
89 flag = 1; # default: do not use the "bias corrected" version
90 else
91 flag = double (flag);
92 if (~ isequal (flag, 0)) && (~ isequal (flag, 1))
93 error ("skewness: FLAG must be 0 or 1");
94 end
56 endif 95 endif
57 96
58 nd = ndims (x); 97 nd = ndims (x);
59 sz = size (x); 98 sz = size (x);
60 if (nargin != 2) 99 if nargin < 3
61 ## Find the first non-singleton dimension. 100 ## Find the first non-singleton dimension.
62 (dim = find (sz > 1, 1)) || (dim = 1); 101 (dim = find (sz > 1, 1)) || (dim = 1);
63 else 102 else
64 if (!(isscalar (dim) && dim == fix (dim)) 103 if (!(isscalar (dim) && dim == fix (dim)) || !(1 <= dim && dim <= nd))
65 || !(1 <= dim && dim <= nd))
66 error ("skewness: DIM must be an integer and a valid dimension"); 104 error ("skewness: DIM must be an integer and a valid dimension");
67 endif 105 endif
68 endif 106 endif
69 107
70 n = sz(dim); 108 n = sz(dim);
71 sz(dim) = 1; 109 sz(dim) = 1;
72 x = center (x, dim); # center also promotes integer to double for next line 110
73 retval = zeros (sz, class (x)); 111 ## In the following sequence of computations, if x is of class single,
74 s = std (x, [], dim); 112 ## then the result y is also of class single
75 idx = find (s > 0); 113 x = center (x, dim);
76 x = sum (x .^ 3, dim); 114 s = std (x, flag, dim); # use the biased estimator of the variance
77 retval(idx) = x(idx) ./ (n * s(idx) .^ 3); 115 y = sum (x .^ 3, dim);
116 y = y ./ (n * s .^ 3);
117 y(s <= 0) = nan;
118
119 ## Apply bias correction to the third central sample moment
120 if flag == 0
121 if n > 2
122 y = y * (n * n ) / ((n - 1) * (n - 2));
123 else
124 y(:) = nan;
125 end
126 endif
78 127
79 endfunction 128 endfunction
80 129
81 130
82 %!assert (skewness ([-1,0,1]), 0) 131 %!assert (skewness ([-1, 0, 1]), 0)
83 %!assert (skewness ([-2,0,1]) < 0) 132 %!assert (skewness ([-2, 0, 1]) < 0)
84 %!assert (skewness ([-1,0,2]) > 0) 133 %!assert (skewness ([-1, 0, 2]) > 0)
85 %!assert (skewness ([-3,0,1]) == -1*skewness ([-1,0,3])) 134 %!assert (skewness ([-3, 0, 1]) == -1 * skewness ([-1, 0, 3]))
135 %!assert (skewness (ones (3, 5)), nan (1, 5))
136
86 %!test 137 %!test
87 %! x = [0; 0; 0; 1]; 138 %! x = [0; 0; 0; 1];
88 %! y = [x, 2*x]; 139 %! y = [x, 2*x];
89 %! assert(all (abs (skewness (y) - [0.75, 0.75]) < sqrt (eps))); 140 %! assert (skewness (y), 1.154700538379251 * [1 1], 1e-13)
90 141
91 %!assert (skewness (single (1)), single (0)) 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], 1, 2), 1.051328089232020 * [1; 1], 1e-15)
144 %!assert (skewness ([1:5 10; 1:5 10], [], 2), 1.051328089232020 * [1; 1], 1e-15)
92 145
93 %% Test input validation 146 ## Test behaviour on single input
147 %! assert (skewness (single ([1:5 10])), single (1.0513283))
148 %! assert (skewness (single ([1 2]), 0), single (nan))
149
150 ## Test input validation
94 %!error skewness () 151 %!error skewness ()
95 %!error skewness (1, 2, 3) 152 %!error skewness (1, 2, 3)
96 %!error skewness (['A'; 'B']) 153 %!error skewness (['A'; 'B'])
97 %!error skewness (1, ones (2,2)) 154 %!error skewness (1, ones (2,2))
98 %!error skewness (1, 1.5) 155 %!error skewness (1, 1.5)
99 %!error skewness (1, 0) 156 %!error skewness (1, [], 0)
100 %!error skewness (1, 3) 157 %!error skewness (1, 3)
101 158