function f = beta_pdf( x, a, b, q, p )
% f = beta_pdf( x, a, b, q, p )
% compute the PDF of the beta probability distribution function
% in terms of it's min and max values (a,b)) and exponents (p,q) 

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

 f = 1/beta(q,p)  * (x-a).^(q-1) .* (b-x).^(p-1) / (b-a)^(q+p-1);

 f(find(x<a)) = 0;
 f(find(x>b)) = 0;

% beta_pdf ---------------------------------------------------- 13 Feb 2024

