function [ X, S ] = mybetarnd( M, N, P, Q, S )
%MYBETARND(M,N,P,Q,S) Draws from a beta distribution using acceptance-
% rejection. This function uses myrand to generate standard uniform
% random numbers and acceptance-rejection based on these uniform random
% numbers to sample from a beta distribution. The parameters M and N
% specify the size of the matrix of (independent) beta random numbers
% which is returned, P and Q are the parameters of the distribution. S is
% the seed which is needed (and changed) by myrand.
% prefactor 1/C with beta function cancelled out
a = ((P+Q-2)^(P+Q-2)) / (((P-1)^(P-1)) * ((Q-1)^(Q-1)));
X = zeros(M,N);
for i=1:M
for j=1:N
[Y, S] = myrand(1, 1, S);
[U, S] = myrand(1, 1, S);
while U > a * (Y^(P-1)) * ((1-Y)^(Q-1))
[Y, S] = myrand(1, 1, S);
[U, S] = myrand(1, 1, S);
end
X(i,j) = Y;
end
end
end