|
7017
|
1 ## Copyright (C) 1995, 1996, 1997, 1998, 2000, 2002, 2003, 2004, 2005, |
|
|
2 ## 2006, 2007 Kurt Hornik |
|
3426
|
3 ## |
|
3922
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
|
7 ## under the terms of the GNU General Public License as published by |
|
7016
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3426
|
10 ## |
|
3922
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
3200
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
3426
|
14 ## General Public License for more details. |
|
|
15 ## |
|
3200
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3200
|
19 |
|
3453
|
20 ## -*- texinfo -*- |
|
4855
|
21 ## @deftypefn {Function File} {} moment (@var{x}, @var{p}, @var{opt}, @var{dim}) |
|
3453
|
22 ## If @var{x} is a vector, compute the @var{p}-th moment of @var{x}. |
|
3200
|
23 ## |
|
3453
|
24 ## If @var{x} is a matrix, return the row vector containing the |
|
|
25 ## @var{p}-th moment of each column. |
|
3200
|
26 ## |
|
|
27 ## With the optional string opt, the kind of moment to be computed can |
|
3453
|
28 ## be specified. If opt contains @code{"c"} or @code{"a"}, central |
|
|
29 ## and/or absolute moments are returned. For example, |
|
|
30 ## |
|
|
31 ## @example |
|
|
32 ## moment (x, 3, "ac") |
|
|
33 ## @end example |
|
|
34 ## |
|
|
35 ## @noindent |
|
|
36 ## computes the third central absolute moment of @var{x}. |
|
4844
|
37 ## |
|
|
38 ## If the optional argument @var{dim} is supplied, work along dimension |
|
|
39 ## @var{dim}. |
|
3453
|
40 ## @end deftypefn |
|
3200
|
41 |
|
|
42 ## Can easily be made to work for continuous distributions (using quad) |
|
|
43 ## as well, but how does the general case work? |
|
3426
|
44 |
|
5428
|
45 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
|
3456
|
46 ## Description: Compute moments |
|
3426
|
47 |
|
4844
|
48 function m = moment (x, p, opt1, opt2) |
|
3426
|
49 |
|
4844
|
50 if ((nargin < 2) || (nargin > 4)) |
|
6046
|
51 print_usage (); |
|
3200
|
52 endif |
|
3426
|
53 |
|
4844
|
54 need_dim = 0; |
|
|
55 |
|
5090
|
56 if (nargin == 2) |
|
4844
|
57 opt = ""; |
|
|
58 need_dim = 1; |
|
5090
|
59 elseif (nargin == 3) |
|
5443
|
60 if (ischar (opt1)) |
|
4844
|
61 opt = opt1; |
|
|
62 need_dim = 1; |
|
|
63 else |
|
|
64 dim = opt1; |
|
|
65 opt = ""; |
|
|
66 endif |
|
5090
|
67 elseif (nargin == 4) |
|
5443
|
68 if (ischar (opt1)) |
|
4844
|
69 opt = opt1; |
|
|
70 dim = opt2; |
|
5443
|
71 elseif (ischar (opt2)) |
|
4844
|
72 opt = opt2; |
|
|
73 dim = opt1; |
|
|
74 else |
|
6046
|
75 error ("moment: expecting opt to be a string"); |
|
4844
|
76 endif |
|
|
77 else |
|
6046
|
78 print_usage (); |
|
3200
|
79 endif |
|
3426
|
80 |
|
4844
|
81 if (need_dim) |
|
|
82 t = find (size (x) != 1); |
|
|
83 if (isempty (t)) |
|
|
84 dim = 1; |
|
|
85 else |
|
|
86 dim = t(1); |
|
|
87 endif |
|
|
88 endif |
|
|
89 |
|
|
90 sz = size (x); |
|
|
91 n = sz (dim); |
|
|
92 |
|
5090
|
93 if (numel (x) < 1) |
|
4844
|
94 error ("moment: x must not be empty"); |
|
3200
|
95 endif |
|
3426
|
96 |
|
5568
|
97 if any (opt == "c") |
|
|
98 rng = ones (1, length (sz)); |
|
|
99 rng(dim) = sz(dim); |
|
|
100 x = x - repmat (sum (x, dim), rng) / n; |
|
|
101 endif |
|
|
102 if any (opt == "a") |
|
|
103 x = abs (x); |
|
|
104 endif |
|
4844
|
105 |
|
5568
|
106 m = sum (x .^ p, dim) / n; |
|
3426
|
107 |
|
3200
|
108 endfunction |
