annotate scripts/statistics/distributions/frnd.m @ 7032:4d4197ffb09d

[project @ 2007-10-15 16:31:55 by jwe]
author jwe
date Mon, 15 Oct 2007 16:31:55 +0000
parents a1dbe9d80eee
children 3306cfcb856e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
1 ## Copyright (C) 1995, 1996, 1997, 2005, 2006, 2007 Kurt Hornik
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
2 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
3 ## This file is part of Octave.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
4 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
6 ## under the terms of the GNU General Public License as published by
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7007
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7007
diff changeset
8 ## your option) any later version.
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
9 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
13 ## General Public License for more details.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
14 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7007
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7007
diff changeset
17 ## <http://www.gnu.org/licenses/>.
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
18
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
19 ## -*- texinfo -*-
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
20 ## @deftypefn {Function File} {} frnd (@var{m}, @var{n}, @var{r}, @var{c})
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
21 ## @deftypefnx {Function File} {} frnd (@var{m}, @var{n}, @var{sz})
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
22 ## Return an @var{r} by @var{c} matrix of random samples from the F
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
23 ## distribution with @var{m} and @var{n} degrees of freedom. Both
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
24 ## @var{m} and @var{n} must be scalar or of size @var{r} by @var{c}.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
25 ## If @var{sz} is a vector the random samples are in a matrix of
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
26 ## size @var{sz}.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
27 ##
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
28 ## If @var{r} and @var{c} are omitted, the size of the result matrix is
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
29 ## the common size of @var{m} and @var{n}.
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
30 ## @end deftypefn
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
31
5428
2a16423e4aa0 [project @ 2005-08-23 18:38:27 by jwe]
jwe
parents: 5411
diff changeset
32 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
33 ## Description: Random deviates from the F distribution
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
34
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
35 function rnd = frnd (m, n, r, c)
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
36
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
37 if (nargin > 1)
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
38 if (!isscalar(m) || !isscalar(n))
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
39 [retval, m, n] = common_size (m, n);
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
40 if (retval > 0)
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
41 error ("frnd: m and n must be of common size or scalar");
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
42 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
43 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
44 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
45
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
46
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
47 if (nargin == 4)
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
48 if (! (isscalar (r) && (r > 0) && (r == round (r))))
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
49 error ("frnd: r must be a positive integer");
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
50 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
51 if (! (isscalar (c) && (c > 0) && (c == round (c))))
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
52 error ("frnd: c must be a positive integer");
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
53 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
54 sz = [r, c];
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
55
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
56 if (any (size (m) != 1) &&
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
57 ((length (size (m)) != length (sz)) || any (size (m) != sz)))
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
58 error ("frnd: m and n must be scalar or of size [r,c]");
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
59 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
60 elseif (nargin == 3)
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
61 if (isscalar (r) && (r > 0))
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
62 sz = [r, r];
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
63 elseif (isvector(r) && all (r > 0))
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
64 sz = r(:)';
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
65 else
7007
6304d9ea0a30 [project @ 2007-10-11 16:26:36 by jwe]
jwe
parents: 6339
diff changeset
66 error ("frnd: r must be a positive integer or vector");
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
67 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
68
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
69 if (any (size (m) != 1) &&
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
70 ((length (size (m)) != length (sz)) || any (size (m) != sz)))
5411
bee21f388110 [project @ 2005-07-13 17:53:44 by jwe]
jwe
parents: 5410
diff changeset
71 error ("frnd: m and n must be scalar or of size sz");
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
72 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
73 elseif (nargin == 2)
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
74 sz = size(a);
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
75 else
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5482
diff changeset
76 print_usage ();
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
77 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
78
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
79
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
80 if (isscalar (m) && isscalar (n))
7032
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
81 if (isinf (m) || isinf (n))
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
82 if (isinf (m))
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
83 rnd = ones (sz);
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
84 else
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
85 rnd = 2 ./ m .* randg(m / 2, sz);
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
86 endif
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
87 if (! isinf (n))
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
88 rnd = 0.5 .* n .* rnd ./ randg (n / 2, sz);
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
89 endif
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
90 elseif ((m > 0) && (m < Inf) && (n > 0) && (n < Inf))
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
91 rnd = n ./ m .* randg (m / 2, sz) ./ randg (n / 2, sz);
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
92 else
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
93 rnd = NaN * ones (sz);
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
94 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
95 else
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
96 rnd = zeros (sz);
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
97
7032
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
98 k = find (isinf(m) | isinf(n));
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
99 if (any (k))
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
100 rnd (k) = 1;
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
101 k2 = find (!isinf(m) & isinf(n));
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
102 rnd (k2) = 2 ./ m(k2) .* randg (m(k2) ./ 2, size(k2));
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
103 k2 = find (isinf(m) & !isinf(n));
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
104 rnd (k2) = 0.5 .* n(k2) .* rnd(k2) ./ randg (n(k2) ./ 2, size(k2));
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
105 endif
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
106
4d4197ffb09d [project @ 2007-10-15 16:31:55 by jwe]
jwe
parents: 7017
diff changeset
107 k = find (!(m > 0) | !(n > 0));
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
108 if (any (k))
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
109 rnd(k) = NaN;
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
110 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
111
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
112 k = find ((m > 0) & (m < Inf) &
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
113 (n > 0) & (n < Inf));
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
114 if (any (k))
6339
20c48710b2c7 [project @ 2007-02-22 18:22:24 by dbateman]
dbateman
parents: 6046
diff changeset
115 rnd(k) = n(k) ./ m(k) .* randg(m(k)./2,size(k)) ./ randg(n(k)./2,size(k));
5410
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
116 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
117 endif
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
118
56e066f5efc1 [project @ 2005-07-13 17:43:35 by jwe]
jwe
parents:
diff changeset
119 endfunction