nsim<-10000;n<-6;p<-.3;
y<-rgamma(nsim,n,p/(1-p));x<-rpois(nsim,y);
hist(x,main="Negative Binomial",freq=F,col="green",breaks=40)
par(new=T)
lines(1:50,dnbinom(1:50,n,p))