function NegativeBinomialDistribution(p,k,sims)
% negative binomial distribution with parameters p and k


Trials = nbinrnd(k,p,1, sims)+k;

n = max(Trials);


for i = 1:n+1
    Dist(i) = (numel(find(Trials==i)))/sims;
end

RR = max(Dist);

plot(1:n+1, Dist, 's', 'MarkerSize', 8)
axis([0 n+1 0 RR+.05])