www.pudn.com > HMM1.zip > beta_sample.m


function r = betarnd(a,b,m,n);
>BETARND Random matrices from beta distribution.
> R = BETARND(A,B) returns a matrix of random numbers chosen
> from the beta distribution with parameters A and B.
> The size of R is the common size of A and B if both are matrices.
> If either parameter is a scalar, the size of R is the size of the other
> parameter. Alternatively, R = BETARND(A,B,M,N) returns an M by N matrix.

> Reference:
> [1] L. Devroye, "Non-Uniform Random Variate Generation",
> Springer-Verlag, 1986

> Copyright (c) 1993-98 by The MathWorks, Inc.
> $Revision: 1.1 $ $Date: 2005/04/26 02:29:17 $

if nargin < 2,
error('Requires at least two input arguments');
end

if nargin == 2
[errorcode rows columns] = rndcheck(2,2,a,b);
end

if nargin == 3
[errorcode rows columns] = rndcheck(3,2,a,b,m);
end

if nargin == 4
[errorcode rows columns] = rndcheck(4,2,a,b,m,n);
end

if errorcode > 0
error('Size information is inconsistent.');
end

r = zeros(rows,columns);

> Use Theorem 4.1, case A (Devroye, page 430) to derive beta
> random numbers as a ratio of gamma random numbers.
if prod(size(a)) == 1
a1 = a(ones(rows,1),ones(columns,1));
g1 = gamrnd(a1,1);
else
g1 = gamrnd(a,1);
end
if prod(size(b)) == 1
b1 = b(ones(rows,1),ones(columns,1));
g2 = gamrnd(b1,1);
else
g2 = gamrnd(b,1);
end
r = g1 ./ (g1 + g2);

> Return NaN if b is not positive.
if any(any(b <= 0));
if prod(size(b) == 1)
tmp = NaN;
r = tmp(ones(rows,columns));
else
k = find(b <= 0);
tmp = NaN;
r(k) = tmp(ones(size(k)));
end
end

> Return NaN if a is not positive.
if any(any(a <= 0));
if prod(size(a) == 1)
tmp = NaN;
r = tmp(ones(rows,columns));
else
k = find(a <= 0);
tmp = NaN;
r(k) = tmp(ones(size(k)));
end
end