function y = nbinom(m, p, k)
% NBINOM y = nbinom(m, p, k) Negative binomial probabilities
% Version of 12/10/92
% Probability the mth success occurs on the kth trial
% m a positive integer;  p a probability
% k a matrix of integers greater than or equal to m
% y = P(X=k) (a matrix of the same dimensions as k)
q = 1 - p;
y = ((p^m)/gamma(m)).*(q.^(k - m)).*gamma(k)./gamma(k - m + 1);