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