function y = binopdf(x,n,p)
% compute the probability mass of the count x 
% occurring as a sum of n binary events, with event probability p
% Inputs:  vector x, 
% n should be either same size as x, or single number
% p should be either same size as x, or single number
% OUTPUTS
% probability mass of the count x.
 nk = gammaln(n + 1) - gammaln(x + 1) - gammaln(n - x + 1);
 lny = nk + x.*log(p)+ (n - x).*log(1 - p);
 y = exp(lny);