alpha = 2;beta = 0.5;n=1000;N=10000;
x = gamrnd(alpha,beta,n,N);
m1 = mean(x);
m2 = mean(x.^2);
alphahat = m1.^2./(m2-m1.^2);
betahat = (m2-m1.^2)./m1;
subplot(1,2,1)
d = 0.1;
b = d/2+[0:d:20];
n = histc(alphahat,b);
bar(b,n/(sum(n)*d),1);
title('alphahat');
subplot(1,2,2)
d = 0.05;
b = d/2+[0:d:3];
n = histc(betahat,b);
bar(b,n/(sum(n)*d),1);
title('betahat');

