r=16;p=0.85;n=10;N=100000;
x = nbinrnd(r,p,n,N);
m1 = mean(x);
m2 = mean(x.^2);
phat = m1./(m2-m1.^2);
rhat = m1.*phat./(1-phat);
subplot(1,2,1)
d = 0.1;
b = d/2+[0:d:2];
n = histc(phat,b);
bar(b,n/(sum(n)*d),1);
title('phat');
subplot(1,2,2)
d = 0.05;
b = d/2+[0:d:30];
n = histc(rhat,b);
bar(b,n/(sum(n)*d),1);
title('rhat');

