n=mean(T)/var(T);
[Y,X]=hist(T,1000);
bar(X,Y)
hold on
npts=sum(Y);
dx=(X(end)-X(1))/(length(X)-1);
Z=(npts*n*dx)*2*(n*X/2).^(n-1).*besselk(0,n*X)/gamma(n/2)^2;
plot(X,Z,'r')

