1######################################################################## 2## 3## Copyright (C) 1995-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {} moment (@var{x}, @var{p}) 28## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{type}) 29## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{dim}) 30## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{type}, @var{dim}) 31## @deftypefnx {} {} moment (@var{x}, @var{p}, @var{dim}, @var{type}) 32## Compute the @var{p}-th central moment of the vector @var{x}. 33## 34## The @var{p}-th central moment of @var{x} is defined as: 35## 36## @tex 37## $$ 38## {\sum_{i=1}^N (x_i - \bar{x})^p \over N} 39## $$ 40## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of elements of @var{x}. 41## 42## 43## @end tex 44## @ifnottex 45## 46## @example 47## @group 48## 1/N SUM_i (@var{x}(i) - mean(@var{x}))^@var{p} 49## @end group 50## @end example 51## 52## @noindent 53## where @math{N} is the length of the @var{x} vector. 54## 55## @end ifnottex 56## 57## If @var{x} is a matrix, return the row vector containing the @var{p}-th 58## central moment of each column. 59## 60## If the optional argument @var{dim} is given, operate along this dimension. 61## 62## The optional string @var{type} specifies the type of moment to be computed. 63## Valid options are: 64## 65## @table @asis 66## @item @qcode{"c"} 67## Central Moment (default). 68## 69## @item @qcode{"a"} 70## @itemx @qcode{"ac"} 71## Absolute Central Moment. The moment about the mean ignoring sign 72## defined as 73## @tex 74## $$ 75## {\sum_{i=1}^N {\left| x_i - \bar{x} \right|}^p \over N} 76## $$ 77## @end tex 78## @ifnottex 79## 80## @example 81## @group 82## 1/N SUM_i (abs (@var{x}(i) - mean(@var{x})))^@var{p} 83## @end group 84## @end example 85## 86## @end ifnottex 87## 88## @item @qcode{"r"} 89## Raw Moment. The moment about zero defined as 90## 91## @tex 92## $$ 93## {\rm moment} (x) = { \sum_{i=1}^N {x_i}^p \over N } 94## $$ 95## @end tex 96## @ifnottex 97## 98## @example 99## @group 100## moment (@var{x}) = 1/N SUM_i @var{x}(i)^@var{p} 101## @end group 102## @end example 103## 104## @end ifnottex 105## 106## @item @nospell{@qcode{"ar"}} 107## Absolute Raw Moment. The moment about zero ignoring sign defined as 108## @tex 109## $$ 110## {\sum_{i=1}^N {\left| x_i \right|}^p \over N} 111## $$ 112## @end tex 113## @ifnottex 114## 115## @example 116## @group 117## 1/N SUM_i ( abs (@var{x}(i)) )^@var{p} 118## @end group 119## @end example 120## 121## @end ifnottex 122## @end table 123## 124## If both @var{type} and @var{dim} are given they may appear in any order. 125## @seealso{var, skewness, kurtosis} 126## @end deftypefn 127 128## Can easily be made to work for continuous distributions (using quad) 129## as well, but how does the general case work? 130 131function m = moment (x, p, opt1, opt2) 132 133 if (nargin < 2 || nargin > 4) 134 print_usage (); 135 endif 136 137 if (! (isnumeric (x) || islogical (x)) || isempty (x)) 138 error ("moment: X must be a non-empty numeric matrix or vector"); 139 endif 140 141 if (! (isnumeric (p) && isscalar (p))) 142 error ("moment: P must be a numeric scalar"); 143 endif 144 145 need_dim = false; 146 147 if (nargin == 2) 148 type = ""; 149 need_dim = true; 150 elseif (nargin == 3) 151 if (ischar (opt1)) 152 type = opt1; 153 need_dim = true; 154 else 155 dim = opt1; 156 type = ""; 157 endif 158 elseif (nargin == 4) 159 if (ischar (opt1)) 160 type = opt1; 161 dim = opt2; 162 elseif (ischar (opt2)) 163 type = opt2; 164 dim = opt1; 165 else 166 error ("moment: TYPE must be a string"); 167 endif 168 endif 169 170 nd = ndims (x); 171 sz = size (x); 172 if (need_dim) 173 ## Find the first non-singleton dimension. 174 (dim = find (sz > 1, 1)) || (dim = 1); 175 else 176 if (! (isscalar (dim) && dim == fix (dim) && dim > 0)) 177 error ("moment: DIM must be an integer and a valid dimension"); 178 endif 179 endif 180 181 n = size (x, dim); 182 183 if (! any (type == "r")) 184 x = center (x, dim); 185 endif 186 if (any (type == "a")) 187 x = abs (x); 188 endif 189 190 m = sum (x .^ p, dim) / n; 191 192endfunction 193 194 195%!shared x 196%! x = rand (10); 197%!assert (moment (x,1), mean (center (x)), eps) 198%!assert (moment (x,2), meansq (center (x)), eps) 199%!assert (moment (x,1,2), mean (center (x, 2), 2), eps) 200%!assert (moment (x,1,"a"), mean (abs (center (x))), eps) 201%!assert (moment (x,1,"r"), mean (x), eps) 202%!assert (moment (x,1,"ar"), mean (abs (x)), eps) 203 204%!assert (moment (single ([1 2 3]), 1, "r"), single (2)) 205 206%!assert (moment (1, 2, 4), 0) 207 208## Test input validation 209%!error moment () 210%!error moment (1) 211%!error moment (1, 2, 3, 4, 5) 212%!error <X must be a non-empty numeric matrix> moment (['A'; 'B'], 2) 213%!error <X must be a non-empty numeric matrix> moment (ones (2,0,3), 2) 214%!error <P must be a numeric scalar> moment (1, true) 215%!error <P must be a numeric scalar> moment (1, ones (2,2)) 216%!error <TYPE must be a string> moment (1, 2, 3, 4) 217%!error <DIM must be an integer and a valid dimension> moment (1, 2, ones (2,2)) 218%!error <DIM must be an integer and a valid dimension> moment (1, 2, 1.5) 219