function mbetaden = fbeta(mx,palpha,pbeta) 
% FBETA, Evaluates Beta(palpha,pbeta) probability density function
% Inputs:
%      mx     - matrix of arguments at which to evaluate density
%      palpha  - scalar parameter alpha, assumed > 0
%      pbeta   - scalar parameter beta, assumed > 0
% Output:
%   mbetaden  - matrix of beta density values, corresponding to
%                       entries of mx
% Note:
%     mean of the Beta(pa,p) distribution is:
%                    E(X) = pa / (pa + pb)
%     variance of the Beta(palpha,pbeta) distribution is:
%       var(X) = (pa * pb) / ((pa + pb)^2 * (pa + pb + 1))
%     the cdf of the Beta(pa,pb) distribution can be
%          calculated using:  BETAINC

%  First get constant multiplier
coeff = 1 / beta(palpha,pbeta) ;

%  Find locations where entries of mx are outside (0,1)
flagout = (mx <= 0)|(mx >= 1) ;
          %  Matrix is 1 where are outside range, 0 elsewhere
          %  Recall | is "or"
flagin = 1 - flagout ;
          %  Matrix is 1 where are inside range, 0 elsewhere

%  Do main calculation
mbetaden = (mx .* flagin + flagout).^(palpha - 1) ;
mbetaden = mbetaden .* ((1 - mx) .* flagin + flagout).^(pbeta - 1) ;
          %  These avoid problems with shifting to complex numbers,
          %      by taking (negative)^(power).
mbetaden = mbetaden .* flagin ;
          %  replace values where f(x) = 0  by  0.
mbetaden = coeff * mbetaden ;
