|
wolffd@0
|
1 function r = betarnd(a,b,m,n);
|
|
wolffd@0
|
2 %BETARND Random matrices from beta distribution.
|
|
wolffd@0
|
3 % R = BETARND(A,B) returns a matrix of random numbers chosen
|
|
wolffd@0
|
4 % from the beta distribution with parameters A and B.
|
|
wolffd@0
|
5 % The size of R is the common size of A and B if both are matrices.
|
|
wolffd@0
|
6 % If either parameter is a scalar, the size of R is the size of the other
|
|
wolffd@0
|
7 % parameter. Alternatively, R = BETARND(A,B,M,N) returns an M by N matrix.
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 % Reference:
|
|
wolffd@0
|
10 % [1] L. Devroye, "Non-Uniform Random Variate Generation",
|
|
wolffd@0
|
11 % Springer-Verlag, 1986
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 % Copyright (c) 1993-98 by The MathWorks, Inc.
|
|
wolffd@0
|
14 % $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 if nargin < 2,
|
|
wolffd@0
|
17 error('Requires at least two input arguments');
|
|
wolffd@0
|
18 end
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 if nargin == 2
|
|
wolffd@0
|
21 [errorcode rows columns] = rndcheck(2,2,a,b);
|
|
wolffd@0
|
22 end
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 if nargin == 3
|
|
wolffd@0
|
25 [errorcode rows columns] = rndcheck(3,2,a,b,m);
|
|
wolffd@0
|
26 end
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28 if nargin == 4
|
|
wolffd@0
|
29 [errorcode rows columns] = rndcheck(4,2,a,b,m,n);
|
|
wolffd@0
|
30 end
|
|
wolffd@0
|
31
|
|
wolffd@0
|
32 if errorcode > 0
|
|
wolffd@0
|
33 error('Size information is inconsistent.');
|
|
wolffd@0
|
34 end
|
|
wolffd@0
|
35
|
|
wolffd@0
|
36 r = zeros(rows,columns);
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 % Use Theorem 4.1, case A (Devroye, page 430) to derive beta
|
|
wolffd@0
|
39 % random numbers as a ratio of gamma random numbers.
|
|
wolffd@0
|
40 if prod(size(a)) == 1
|
|
wolffd@0
|
41 a1 = a(ones(rows,1),ones(columns,1));
|
|
wolffd@0
|
42 g1 = gamrnd(a1,1);
|
|
wolffd@0
|
43 else
|
|
wolffd@0
|
44 g1 = gamrnd(a,1);
|
|
wolffd@0
|
45 end
|
|
wolffd@0
|
46 if prod(size(b)) == 1
|
|
wolffd@0
|
47 b1 = b(ones(rows,1),ones(columns,1));
|
|
wolffd@0
|
48 g2 = gamrnd(b1,1);
|
|
wolffd@0
|
49 else
|
|
wolffd@0
|
50 g2 = gamrnd(b,1);
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52 r = g1 ./ (g1 + g2);
|
|
wolffd@0
|
53
|
|
wolffd@0
|
54 % Return NaN if b is not positive.
|
|
wolffd@0
|
55 if any(any(b <= 0));
|
|
wolffd@0
|
56 if prod(size(b) == 1)
|
|
wolffd@0
|
57 tmp = NaN;
|
|
wolffd@0
|
58 r = tmp(ones(rows,columns));
|
|
wolffd@0
|
59 else
|
|
wolffd@0
|
60 k = find(b <= 0);
|
|
wolffd@0
|
61 tmp = NaN;
|
|
wolffd@0
|
62 r(k) = tmp(ones(size(k)));
|
|
wolffd@0
|
63 end
|
|
wolffd@0
|
64 end
|
|
wolffd@0
|
65
|
|
wolffd@0
|
66 % Return NaN if a is not positive.
|
|
wolffd@0
|
67 if any(any(a <= 0));
|
|
wolffd@0
|
68 if prod(size(a) == 1)
|
|
wolffd@0
|
69 tmp = NaN;
|
|
wolffd@0
|
70 r = tmp(ones(rows,columns));
|
|
wolffd@0
|
71 else
|
|
wolffd@0
|
72 k = find(a <= 0);
|
|
wolffd@0
|
73 tmp = NaN;
|
|
wolffd@0
|
74 r(k) = tmp(ones(size(k)));
|
|
wolffd@0
|
75 end
|
|
wolffd@0
|
76 end
|
