|
3200
|
1 ## Copyright (C) 1995, 1996, 1997 Kurt Hornik |
|
3426
|
2 ## |
|
3922
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
3200
|
7 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
8 ## any later version. |
|
3426
|
9 ## |
|
3922
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
3200
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
3426
|
13 ## General Public License for more details. |
|
|
14 ## |
|
3200
|
15 ## You should have received a copy of the GNU General Public License |
|
3922
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA |
|
|
18 ## 02111-1307, USA. |
|
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 |
|
3456
|
45 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> |
|
|
46 ## Description: Compute moments |
|
3426
|
47 |
|
4844
|
48 function m = moment (x, p, opt1, opt2) |
|
3426
|
49 |
|
4844
|
50 if ((nargin < 2) || (nargin > 4)) |
|
|
51 usage ("moment (x, p, type, dim)"); |
|
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) |
|
4844
|
60 if (isstr (opt1)) |
|
|
61 opt = opt1; |
|
|
62 need_dim = 1; |
|
|
63 else |
|
|
64 dim = opt1; |
|
|
65 opt = ""; |
|
|
66 endif |
|
5090
|
67 elseif (nargin == 4) |
|
4844
|
68 if (isstr (opt1)) |
|
|
69 opt = opt1; |
|
|
70 dim = opt2; |
|
|
71 elseif (isstr (opt2)) |
|
|
72 opt = opt2; |
|
|
73 dim = opt1; |
|
|
74 else |
|
|
75 usage ("moment: expecting opt to be a string"); |
|
|
76 endif |
|
|
77 else |
|
|
78 usage ("moment (x, p, dim, opt) or moment (x, p, dim, opt)"); |
|
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 |
|
4844
|
97 tmp = warn_str_to_num; |
|
|
98 unwind_protect |
|
|
99 warn_str_to_num = 0; |
|
|
100 if any (opt == "c") |
|
|
101 rng = ones(1, length (sz)); |
|
|
102 rng (dim) = sz (dim); |
|
|
103 x = x - repmat (sum (x, dim), rng) / n; |
|
|
104 endif |
|
|
105 if any (opt == "a") |
|
|
106 x = abs (x); |
|
|
107 endif |
|
|
108 unwind_protect_cleanup |
|
|
109 warn_str_to_num = tmp; |
|
|
110 end_unwind_protect |
|
|
111 |
|
|
112 m = sum(x .^ p, dim) / n; |
|
3426
|
113 |
|
3200
|
114 endfunction |
