## -*- Mode: octave -*-
function p=betapdf(x, alpha, beta)
  z = lgamma(alpha + beta) - lgamma(alpha) - lgamma(beta);
  ya = log(x) .* (alpha - 1);
  ya(x==0) = 0;
  yb = log(1 - x) .* (beta - 1);
  p = exp(z + ya + yb);
  p(x==0) = 0;
  p(x==1) = 0;
end
