|
wolffd@0
|
1 function r = gamrnd(a,b,m,n);
|
|
wolffd@0
|
2 %GAMRND Random matrices from gamma distribution.
|
|
wolffd@0
|
3 % R = GAMRND(A,B) returns a matrix of random numbers chosen
|
|
wolffd@0
|
4 % from the gamma 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 = GAMRND(A,B,M,N) returns an M by N matrix.
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % Some references refer to the gamma distribution
|
|
wolffd@0
|
10 % with a single parameter. This corresponds to GAMRND
|
|
wolffd@0
|
11 % with B = 1. (See Devroye, pages 401-402.)
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 % GAMRND uses a rejection or an inversion method depending on the
|
|
wolffd@0
|
14 % value of A.
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 % References:
|
|
wolffd@0
|
17 % [1] L. Devroye, "Non-Uniform Random Variate Generation",
|
|
wolffd@0
|
18 % Springer-Verlag, 1986
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 % B.A. Jones 2-1-93
|
|
wolffd@0
|
21 % Copyright (c) 1993-98 by The MathWorks, Inc.
|
|
wolffd@0
|
22 % $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 if nargin < 2,
|
|
wolffd@0
|
25 error('Requires at least two input arguments.');
|
|
wolffd@0
|
26 end
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 if nargin == 2
|
|
wolffd@0
|
30 [errorcode rows columns] = rndcheck(2,2,a,b);
|
|
wolffd@0
|
31 end
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33 if nargin == 3
|
|
wolffd@0
|
34 [errorcode rows columns] = rndcheck(3,2,a,b,m);
|
|
wolffd@0
|
35 end
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37 if nargin == 4
|
|
wolffd@0
|
38 [errorcode rows columns] = rndcheck(4,2,a,b,m,n);
|
|
wolffd@0
|
39 end
|
|
wolffd@0
|
40
|
|
wolffd@0
|
41 if errorcode > 0
|
|
wolffd@0
|
42 error('Size information is inconsistent.');
|
|
wolffd@0
|
43 end
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 % Initialize r to zero.
|
|
wolffd@0
|
46 lth = rows*columns;
|
|
wolffd@0
|
47 r = zeros(lth,1);
|
|
wolffd@0
|
48 a = a(:); b = b(:);
|
|
wolffd@0
|
49
|
|
wolffd@0
|
50 scalara = (length(a) == 1);
|
|
wolffd@0
|
51 if scalara
|
|
wolffd@0
|
52 a = a*ones(lth,1);
|
|
wolffd@0
|
53 end
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 scalarb = (length(b) == 1);
|
|
wolffd@0
|
56 if scalarb
|
|
wolffd@0
|
57 b = b*ones(lth,1);
|
|
wolffd@0
|
58 end
|
|
wolffd@0
|
59
|
|
wolffd@0
|
60 % If a == 1, then gamma is exponential. (Devroye, page 405).
|
|
wolffd@0
|
61 k = find(a == 1);
|
|
wolffd@0
|
62 if any(k)
|
|
wolffd@0
|
63 r(k) = -b(k) .* log(rand(size(k)));
|
|
wolffd@0
|
64 end
|
|
wolffd@0
|
65
|
|
wolffd@0
|
66
|
|
wolffd@0
|
67 k = find(a < 1 & a > 0);
|
|
wolffd@0
|
68 % (Devroye, page 418 Johnk's generator)
|
|
wolffd@0
|
69 if any(k)
|
|
wolffd@0
|
70 c = zeros(lth,1);
|
|
wolffd@0
|
71 d = zeros(lth,1);
|
|
wolffd@0
|
72 c(k) = 1 ./ a(k);
|
|
wolffd@0
|
73 d(k) = 1 ./ (1 - a(k));
|
|
wolffd@0
|
74 accept = k;
|
|
wolffd@0
|
75 while ~isempty(accept)
|
|
wolffd@0
|
76 u = rand(size(accept));
|
|
wolffd@0
|
77 v = rand(size(accept));
|
|
wolffd@0
|
78 x = u .^ c(accept);
|
|
wolffd@0
|
79 y = v .^ d(accept);
|
|
wolffd@0
|
80 k1 = find((x + y) <= 1);
|
|
wolffd@0
|
81 if ~isempty(k1)
|
|
wolffd@0
|
82 e = -log(rand(size(k1)));
|
|
wolffd@0
|
83 r(accept(k1)) = e .* x(k1) ./ (x(k1) + y(k1));
|
|
wolffd@0
|
84 accept(k1) = [];
|
|
wolffd@0
|
85 end
|
|
wolffd@0
|
86 end
|
|
wolffd@0
|
87 r(k) = r(k) .* b(k);
|
|
wolffd@0
|
88 end
|
|
wolffd@0
|
89
|
|
wolffd@0
|
90 % Use a rejection method for a > 1.
|
|
wolffd@0
|
91 k = find(a > 1);
|
|
wolffd@0
|
92 % (Devroye, page 410 Best's algorithm)
|
|
wolffd@0
|
93 bb = zeros(size(a));
|
|
wolffd@0
|
94 c = bb;
|
|
wolffd@0
|
95 if any(k)
|
|
wolffd@0
|
96 bb(k) = a(k) - 1;
|
|
wolffd@0
|
97 c(k) = 3 * a(k) - 3/4;
|
|
wolffd@0
|
98 accept = k;
|
|
wolffd@0
|
99 count = 1;
|
|
wolffd@0
|
100 while ~isempty(accept)
|
|
wolffd@0
|
101 m = length(accept);
|
|
wolffd@0
|
102 u = rand(m,1);
|
|
wolffd@0
|
103 v = rand(m,1);
|
|
wolffd@0
|
104 w = u .* (1 - u);
|
|
wolffd@0
|
105 y = sqrt(c(accept) ./ w) .* (u - 0.5);
|
|
wolffd@0
|
106 x = bb(accept) + y;
|
|
wolffd@0
|
107 k1 = find(x >= 0);
|
|
wolffd@0
|
108 if ~isempty(k1)
|
|
wolffd@0
|
109 z = 64 * (w .^ 3) .* (v .^ 2);
|
|
wolffd@0
|
110 k2 = (z(k1) <= (1 - 2 * (y(k1) .^2) ./ x(k1)));
|
|
wolffd@0
|
111 k3 = k1(find(k2));
|
|
wolffd@0
|
112 r(accept(k3)) = x(k3);
|
|
wolffd@0
|
113 k4 = k1(find(~k2));
|
|
wolffd@0
|
114 k5 = k4(find(log(z(k4)) <= (2*(bb(accept(k4)).*log(x(k4)./bb(accept(k4)))-y(k4)))));
|
|
wolffd@0
|
115 r(accept(k5)) = x(k5);
|
|
wolffd@0
|
116 omit = [k3; k5];
|
|
wolffd@0
|
117 accept(omit) = [];
|
|
wolffd@0
|
118 end
|
|
wolffd@0
|
119 end
|
|
wolffd@0
|
120 r(k) = r(k) .* b(k);
|
|
wolffd@0
|
121 end
|
|
wolffd@0
|
122
|
|
wolffd@0
|
123 % Return NaN if a or b is not positive.
|
|
wolffd@0
|
124 r(b <= 0 | a <= 0) = NaN;
|
|
wolffd@0
|
125
|
|
wolffd@0
|
126 r = reshape(r,rows,columns);
|
