function x = beta_rnd( a, b, q, p, M, N )
% f = beta_rnd( a, b, q, p, M, N )
% compute a sample with a beta probability distribution 
% in terms of it's min and max values (a,b)) and exponents (q,p) 

 if (b < a)
   printf('error: beta_rnd: a = %f, b = %f, a must be less than b \n', a, b );
 end

 X  =  randg(q,M,N);
 Y  =  randg(p,M,N);

 z = X ./ (X + Y);

 x = a + (b-a)*z;

% beta_rnd ---------------------------------------------------- 13 Feb 2024

